Analyses
Nicolas O. Rode
22/02/2021
Main conclusions: Statistical power higher for emergence rate than for number of adults Emergence rate is density dependent with optimum around 100eggs per tube Variation in the number of eggs between 50 and 100 has more effect in emeergence rates than variation between 150 and 250 The SA model should include the number of eggs as a covariate or a factor wit an egg score to increase power
1 Load packages
library(lattice)
library(lme4)## Loading required package: Matrix2 Performance
2.1 Load data
data <- read.table(file=here::here("data", "DATACOMPLET_PERF.csv"), header=TRUE, sep=",")
head(data)## Generation Experiment Original_environment Population Date_P Date_C_O
## 1 G0 Plasticity WT3 WT3 18/09/2018 19/09/2018
## 2 G0 Plasticity WT3 WT3 18/09/2018 19/09/2018
## 3 G0 Plasticity WT3 WT3 18/09/2018 19/09/2018
## 4 G0 Plasticity WT3 WT3 18/09/2018 19/09/2018
## 5 G0 Plasticity WT3 WT3 18/09/2018 19/09/2018
## 6 G0 Plasticity Blackberry Blackberry45 3/10/18 4/10/18
## Date_C_A Row Column Rack Test_environment Nb_eggs Obs_O Nb_adults Obs_A
## 1 4/10/18 L1 C1 1 Strawberry 41 LO 2 LO
## 2 4/10/18 L1 C2 1 Blackberry 85 LO 32 LO
## 3 4/10/18 L1 C3 1 GF 63 LO 14 LO
## 4 4/10/18 L1 C4 1 Cherry 77 LO 28 LO
## 5 4/10/18 L1 C5 1 Grape 38 LO 19 LO
## 6 19/10/2018 L1 C1 1 Blackberry 0 LO 0 LOstr(data)## 'data.frame': 2150 obs. of 15 variables:
## $ Generation : Factor w/ 2 levels "G0","G2": 1 1 1 1 1 1 1 1 1 1 ...
## $ Experiment : Factor w/ 2 levels "Adaptation","Plasticity": 2 2 2 2 2 2 2 2 2 2 ...
## $ Original_environment: Factor w/ 4 levels "Blackberry","Cherry",..: 4 4 4 4 4 1 1 1 1 1 ...
## $ Population : Factor w/ 26 levels "Blackberry31",..: 26 26 26 26 26 13 13 13 13 13 ...
## $ Date_P : Factor w/ 33 levels "1/10/18","10/7/18",..: 11 11 11 11 11 26 26 26 26 26 ...
## $ Date_C_O : Factor w/ 33 levels "1/8/18","11/7/18",..: 11 11 11 11 11 27 27 27 27 27 ...
## $ Date_C_A : Factor w/ 33 levels "1/12/18","11/10/18",..: 28 28 28 28 28 16 16 16 16 16 ...
## $ Row : Factor w/ 10 levels "L1","L10","L2",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Column : Factor w/ 10 levels "C1","C10","C2",..: 1 3 4 5 6 1 3 4 5 6 ...
## $ Rack : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Test_environment : Factor w/ 5 levels "Blackberry","Cherry",..: 5 1 3 2 4 1 3 2 4 5 ...
## $ Nb_eggs : int 41 85 63 77 38 0 0 3 1 0 ...
## $ Obs_O : Factor w/ 6 levels "","AA","LO","NR",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ Nb_adults : int 2 32 14 28 19 0 0 2 2 0 ...
## $ Obs_A : Factor w/ 6 levels "","AA","CD","JL",..: 5 5 5 5 5 5 5 5 5 5 ...data$Generation## [1] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [25] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [49] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [73] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [97] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [121] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [145] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [169] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [193] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [217] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [241] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [265] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [289] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [313] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [337] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [361] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [385] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [409] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [433] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [457] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [481] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [505] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [529] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [553] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [577] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [601] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [625] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [649] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [673] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [697] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [721] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [745] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [769] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [793] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [817] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [841] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [865] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [889] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [913] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [937] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [961] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [985] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [1009] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [1033] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0
## [1057] G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G0 G2 G2 G2 G2 G2 G2 G2
## [1081] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1105] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1129] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1153] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1177] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1201] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1225] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1249] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1273] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1297] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1321] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1345] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1369] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1393] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1417] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1441] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1465] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1489] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1513] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1537] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1561] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1585] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1609] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1633] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1657] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1681] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1705] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1729] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1753] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1777] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1801] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1825] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1849] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1873] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1897] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1921] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1945] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1969] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [1993] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2017] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2041] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2065] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2089] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2113] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## [2137] G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2 G2
## Levels: G0 G2data <- data[data$Generation=="G2",]
#Remove GF and Grape
data <- data[data$Test_environment!="GF",]
data <- data[data$Test_environment!="Grape",]
data$Test_environment <- factor(data$Test_environment)
data$EggScore <- as.factor(ifelse(data$Nb_eggs<51, 1, ifelse(data$Nb_eggs<101, 2, ifelse(data$Nb_eggs<151, 3, 4))))
data$EggScoreFive <- as.factor(ifelse(data$Nb_eggs<51, 1, ifelse(data$Nb_eggs<101, 2, ifelse(data$Nb_eggs<151, 3, ifelse(data$Nb_eggs<201, 4, 5)))))
data$EggScoreSmall <- as.factor(ifelse(data$Nb_eggs<51, 1, ifelse(data$Nb_eggs<76, 2, ifelse(data$Nb_eggs<101, 3, ifelse(data$Nb_eggs<126, 4, ifelse(data$Nb_eggs<151, 5, 6))))))
## Remove WT3
data <- data[data$Population!="WT3",]
data$Original_environment <- factor(data$Original_environment)
## Remove Cherry52 which has a very high mean
#data <- data[data$Population!="Cherry52",]
data$Population <- factor(data$Population)
data$SA <- as.factor(ifelse(data$Original_environment==data$Test_environment, 1, 0))
data$Emergence_Rate <- data$Nb_adults/data$Nb_eggs
head(data)## Generation Experiment Original_environment Population Date_P
## 1074 G2 Adaptation Cherry Cherry3 27/08/2018
## 1075 G2 Adaptation Cherry Cherry3 27/08/2018
## 1078 G2 Adaptation Cherry Cherry3 27/08/2018
## 1079 G2 Adaptation Cherry Cherry3 27/08/2018
## 1080 G2 Adaptation Cherry Cherry3 27/08/2018
## 1083 G2 Adaptation Cherry Cherry3 27/08/2018
## Date_C_O Date_C_A Row Column Rack Test_environment Nb_eggs Obs_O
## 1074 28/08/2018 12/9/18 L1 C1 1 Cherry 180 LO
## 1075 28/08/2018 12/9/18 L1 C2 1 Blackberry 85 LO
## 1078 28/08/2018 12/9/18 L1 C5 1 Strawberry 67 LO
## 1079 28/08/2018 12/9/18 L2 C1 1 Cherry 106 LO
## 1080 28/08/2018 12/9/18 L2 C2 1 Blackberry 113 LO
## 1083 28/08/2018 12/9/18 L2 C5 1 Strawberry 78 LO
## Nb_adults Obs_A EggScore EggScoreFive EggScoreSmall SA Emergence_Rate
## 1074 7 LO 4 4 6 1 0.03888889
## 1075 3 LO 2 2 3 0 0.03529412
## 1078 0 LO 2 2 2 0 0.00000000
## 1079 2 LO 3 3 4 1 0.01886792
## 1080 0 LO 3 3 4 0 0.00000000
## 1083 0 LO 2 2 3 0 0.00000000## Compute emergence rate
data$Emergence_Rate[data$Emergence_Rate>1] <- rep(1, sum(data$Emergence_Rate>1))2.2 Number of eggs
2.2.1 Explore data
histogram(~Nb_eggs|Original_environment*Test_environment, data=data)
tapply(data$Nb_eggs, list(data$Original_environment, data$Test_environment), mean)## Blackberry Cherry Strawberry
## Blackberry 103.9720 121.8491 98.71963
## Cherry 131.7174 133.8400 103.51064
## Strawberry 119.1739 119.4894 111.565222.2.2 Analysis
m1 <- lm(log(Nb_eggs+1) ~ Population + Test_environment + SA + Original_environment:Test_environment, data=data[data$Generation=="G2",])
summary(m1)##
## Call:
## lm(formula = log(Nb_eggs + 1) ~ Population + Test_environment +
## SA + Original_environment:Test_environment, data = data[data$Generation ==
## "G2", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2545 -0.1778 0.0195 0.2192 1.0845
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error
## (Intercept) 4.5806543 0.1382464
## PopulationBlackberry32 -0.6895997 0.1067135
## PopulationBlackberry33 0.0419347 0.1042064
## PopulationBlackberry34 -0.2084134 0.1524978
## PopulationBlackberry35 -0.1012098 0.1033412
## PopulationBlackberry36 -0.0001252 0.1319869
## PopulationBlackberry37 -0.0843005 0.1106111
## PopulationBlackberry38 0.1732260 0.2076342
## PopulationBlackberry39 0.1099952 0.1195802
## PopulationBlackberry40 -0.2846811 0.1130724
## PopulationBlackberry43 -0.3354278 0.1578266
## PopulationBlackberry44 0.1972677 0.1130724
## PopulationBlackberry45 -0.0814644 0.2285211
## PopulationCherry103 -0.1405432 0.1419273
## PopulationCherry104 -0.8139547 0.1482000
## PopulationCherry3 -0.1536887 0.1750231
## PopulationCherry47 0.1054940 0.1228905
## PopulationCherry50 0.1919031 0.1635145
## PopulationCherry51 0.0426110 0.4444066
## PopulationCherry52 -0.1050116 0.2381451
## PopulationCherry6 0.3390668 0.1889877
## PopulationCherry7 -0.1152803 0.1241927
## PopulationStrawberry42 0.1155414 0.1196040
## PopulationStrawberry44 0.1533210 0.1285429
## PopulationStrawberry53 -0.4392990 0.1155464
## Test_environmentCherry 0.2938095 0.1345022
## Test_environmentStrawberry 0.0549272 0.1060046
## SA1 0.0565484 0.1059980
## Test_environmentBlackberry:Original_environmentCherry 0.3183771 0.1384989
## Test_environmentCherry:Original_environmentCherry 0.0011227 0.1703268
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.2444723 0.1834154
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## t value Pr(>|t|)
## (Intercept) 33.134 < 2e-16 ***
## PopulationBlackberry32 -6.462 2.21e-10 ***
## PopulationBlackberry33 0.402 0.687526
## PopulationBlackberry34 -1.367 0.172268
## PopulationBlackberry35 -0.979 0.327809
## PopulationBlackberry36 -0.001 0.999244
## PopulationBlackberry37 -0.762 0.446294
## PopulationBlackberry38 0.834 0.404470
## PopulationBlackberry39 0.920 0.358043
## PopulationBlackberry40 -2.518 0.012085 *
## PopulationBlackberry43 -2.125 0.033991 *
## PopulationBlackberry44 1.745 0.081590 .
## PopulationBlackberry45 -0.356 0.721609
## PopulationCherry103 -0.990 0.322473
## PopulationCherry104 -5.492 5.98e-08 ***
## PopulationCherry3 -0.878 0.380256
## PopulationCherry47 0.858 0.391010
## PopulationCherry50 1.174 0.241038
## PopulationCherry51 0.096 0.923647
## PopulationCherry52 -0.441 0.659411
## PopulationCherry6 1.794 0.073323 .
## PopulationCherry7 -0.928 0.353677
## PopulationStrawberry42 0.966 0.334436
## PopulationStrawberry44 1.193 0.233458
## PopulationStrawberry53 -3.802 0.000159 ***
## Test_environmentCherry 2.184 0.029337 *
## Test_environmentStrawberry 0.518 0.604548
## SA1 0.533 0.593905
## Test_environmentBlackberry:Original_environmentCherry 2.299 0.021878 *
## Test_environmentCherry:Original_environmentCherry 0.007 0.994743
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 1.333 0.183100
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4262 on 571 degrees of freedom
## Multiple R-squared: 0.3375, Adjusted R-squared: 0.3027
## F-statistic: 9.697 on 30 and 571 DF, p-value: < 2.2e-16anova(m1)## Analysis of Variance Table
##
## Response: log(Nb_eggs + 1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Population 24 45.525 1.89686 10.4403 < 2.2e-16 ***
## Test_environment 2 5.834 2.91691 16.0546 1.647e-07 ***
## SA 1 0.404 0.40421 2.2248 0.1364
## Test_environment:Original_environment 3 1.093 0.36447 2.0060 0.1120
## Residuals 571 103.743 0.18169
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## F test for SA
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))## [1] 1.10903(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.3696212## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))## [1] 0.026534192.3 Number of adults
2.3.1 Explore data with blue
m0 <- lm(Nb_adults~Original_environment*Test_environment, data=data[data$Generation=="G2",])
summary(m0)##
## Call:
## lm(formula = Nb_adults ~ Original_environment * Test_environment,
## data = data[data$Generation == "G2", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.318 -11.287 -2.906 8.834 85.087
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 27.3178 1.5742
## Original_environmentCherry -4.4047 2.8710
## Original_environmentStrawberry -1.6656 2.8710
## Test_environmentCherry -4.1291 2.2315
## Test_environmentStrawberry -14.1589 2.2263
## Original_environmentCherry:Test_environmentCherry 4.5360 4.0059
## Original_environmentStrawberry:Test_environmentCherry 0.8386 4.0479
## Original_environmentCherry:Test_environmentStrawberry 2.7777 4.0450
## Original_environmentStrawberry:Test_environmentStrawberry 4.5502 4.0602
## t value Pr(>|t|)
## (Intercept) 17.353 < 2e-16 ***
## Original_environmentCherry -1.534 0.1255
## Original_environmentStrawberry -0.580 0.5620
## Test_environmentCherry -1.850 0.0648 .
## Test_environmentStrawberry -6.360 4.03e-10 ***
## Original_environmentCherry:Test_environmentCherry 1.132 0.2579
## Original_environmentStrawberry:Test_environmentCherry 0.207 0.8359
## Original_environmentCherry:Test_environmentStrawberry 0.687 0.4925
## Original_environmentStrawberry:Test_environmentStrawberry 1.121 0.2629
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 16.28 on 593 degrees of freedom
## Multiple R-squared: 0.1037, Adjusted R-squared: 0.09163
## F-statistic: 8.578 on 8 and 593 DF, p-value: 4.287e-11tapply(data$Nb_adults[data$Generation=="G2"], list(data$Original_environment[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)## Blackberry Cherry Strawberry
## Blackberry 27.31776 23.18868 13.15888
## Cherry 22.91304 23.32000 11.53191
## Strawberry 25.65217 22.36170 16.04348tapply(data$Nb_adults[data$Generation=="G2"], list(data$Original_environment[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), var)## Blackberry Cherry Strawberry
## Blackberry 402.3698 205.4688 143.9085
## Cherry 506.6589 412.7935 174.1240
## Strawberry 337.0763 133.4968 120.3092range(data$Nb_adults)## [1] 0 108## Check number of eggs and adults
tapply(data$Nb_eggs[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)## Blackberry Cherry Strawberry
## Blackberry31 106.44444 130.33333 118.88889
## Blackberry32 61.23077 84.46154 74.00000
## Blackberry33 125.13333 141.28571 95.66667
## Blackberry34 105.75000 95.66667 83.50000
## Blackberry35 93.60000 123.73333 96.31250
## Blackberry36 118.00000 144.00000 92.83333
## Blackberry37 108.09091 114.45455 101.27273
## Blackberry38 125.50000 141.00000 143.00000
## Blackberry39 133.75000 146.00000 117.12500
## Blackberry40 88.30000 92.10000 80.40000
## Blackberry43 69.33333 90.00000 86.33333
## Blackberry44 129.20000 155.70000 129.70000
## Blackberry45 92.00000 113.00000 110.00000
## Cherry103 129.66667 118.00000 89.83333
## Cherry104 59.80000 61.83333 55.20000
## Cherry3 109.66667 143.66667 78.00000
## Cherry47 142.00000 162.71429 128.23077
## Cherry50 148.33333 185.00000 126.00000
## Cherry51 139.00000 NA NA
## Cherry52 129.50000 96.00000 103.00000
## Cherry6 226.50000 203.50000 124.00000
## Cherry7 137.75000 118.69231 97.50000
## Strawberry42 140.46667 131.81250 146.66667
## Strawberry44 150.60000 142.30000 139.70000
## Strawberry53 89.00000 99.23810 73.09524tapply(data$Nb_adults[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)## Blackberry Cherry Strawberry
## Blackberry31 38.444444 36.777778 19.2222222
## Blackberry32 1.923077 6.307692 0.6923077
## Blackberry33 30.800000 34.000000 16.6666667
## Blackberry34 13.500000 17.666667 8.0000000
## Blackberry35 14.533333 23.533333 5.2500000
## Blackberry36 28.800000 23.166667 14.0000000
## Blackberry37 34.727273 24.090909 8.2727273
## Blackberry38 15.500000 7.000000 1.0000000
## Blackberry39 28.375000 26.750000 19.0000000
## Blackberry40 38.700000 17.300000 23.9000000
## Blackberry43 45.000000 13.750000 11.6666667
## Blackberry44 48.900000 27.300000 25.6000000
## Blackberry45 11.500000 30.000000 2.0000000
## Cherry103 35.666667 41.857143 20.8333333
## Cherry104 8.200000 10.833333 3.0000000
## Cherry3 1.000000 3.000000 0.0000000
## Cherry47 15.833333 18.571429 5.9230769
## Cherry50 63.000000 67.750000 46.0000000
## Cherry51 78.000000 NA NA
## Cherry52 65.500000 39.000000 18.0000000
## Cherry6 5.500000 10.000000 1.3333333
## Cherry7 16.416667 16.076923 9.9166667
## Strawberry42 45.733333 35.000000 26.3333333
## Strawberry44 13.100000 13.400000 9.9000000
## Strawberry53 17.285714 17.000000 11.6190476## Check for the presence of negative correlations
m1 <- lm(log(Nb_adults+1) ~ Population + Test_environment, data=data[data$Generation=="G2",])
data$resid <- residuals(m1)
meanbypopbytestenv <- as.data.frame(tapply(data$resid[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean))
## Cherry ~ Blackberry
plot(meanbypopbytestenv$Blackberry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Cherry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Strawberry ~ Cherry
plot(meanbypopbytestenv$Cherry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Cherry", ylab="Strawberry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Strawberry ~ Blackberry
plot(meanbypopbytestenv$Blackberry~ meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Strawberry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
2.3.2 Explore data with BLUPS
m0 <- lmer(log(Nb_adults+1)~ Population + Test_environment + (0 + Test_environment | Population), data = data[data$Generation=="G2",])## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 1 negative eigenvalues## Correlation between fixed effects
plot(coef(m1), fixef(m0))
summary(m0)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## log(Nb_adults + 1) ~ Population + Test_environment + (0 + Test_environment |
## Population)
## Data: data[data$Generation == "G2", ]
##
## REML criterion at convergence: 1264.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8381 -0.4992 0.1024 0.5908 2.6516
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Population Test_environmentBlackberry 9.254e-07 0.000962
## Test_environmentCherry 8.208e-02 0.286500 0.59
## Test_environmentStrawberry 4.287e-02 0.207058 0.44 -0.13
## Residual 4.404e-01 0.663613
## Number of obs: 602, groups: Population, 25
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 3.55606 0.16116 22.066
## PopulationBlackberry32 -2.50179 0.20907 -11.966
## PopulationBlackberry33 -0.17462 0.20469 -0.853
## PopulationBlackberry34 -1.17436 0.27159 -4.324
## PopulationBlackberry35 -1.04239 0.20424 -5.104
## PopulationBlackberry36 -0.32019 0.24803 -1.291
## PopulationBlackberry37 -0.52201 0.21491 -2.429
## PopulationBlackberry38 -1.36372 0.35317 -3.861
## PopulationBlackberry39 -0.34622 0.22792 -1.519
## PopulationBlackberry40 -0.09950 0.21853 -0.455
## PopulationBlackberry43 -0.48439 0.28519 -1.698
## PopulationBlackberry44 0.16978 0.21853 0.777
## PopulationBlackberry45 -1.09654 0.37593 -2.917
## PopulationCherry103 0.08807 0.24101 0.365
## PopulationCherry104 -1.51454 0.25101 -6.034
## PopulationCherry3 -2.86774 0.29040 -9.875
## PopulationCherry47 -1.07034 0.21085 -5.076
## PopulationCherry50 0.68579 0.27800 2.467
## PopulationCherry51 0.81338 0.68290 1.191
## PopulationCherry52 0.37570 0.37593 0.999
## PopulationCherry6 -1.98132 0.31557 -6.279
## PopulationCherry7 -0.82275 0.21155 -3.889
## PopulationStrawberry42 0.23210 0.20440 1.136
## PopulationStrawberry44 -0.93199 0.21853 -4.265
## PopulationStrawberry53 -0.75587 0.19532 -3.870
## Test_environmentCherry 0.04809 0.09288 0.518
## Test_environmentStrawberry -0.70086 0.08215 -8.531##
## Correlation matrix not shown by default, as p = 27 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 1 negative eigenvaluesBLUP <- ranef(m0)$Population
pairs(BLUP)
## Negative correlation between blups and blues, very weird, no?
plot(BLUP$Test_environmentBlackberry, meanbypopbytestenv$Blackberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentCherry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentStrawberry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
2.3.3 Analysis
## Null model
m0 <- lm(log(Nb_adults+1) ~ Population + Test_environment + Original_environment:Test_environment, data=data[data$Generation=="G2",])
## SA model
m1 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment, data=data[data$Generation=="G2",])
summary(m1)##
## Call:
## lm(formula = log(Nb_adults + 1) ~ Population + Test_environment +
## SA + Original_environment:Test_environment, data = data[data$Generation ==
## "G2", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.50160 -0.35356 0.05873 0.42543 1.98747
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error
## (Intercept) 3.20150 0.22142
## PopulationBlackberry32 -2.35665 0.17091
## PopulationBlackberry33 -0.12012 0.16690
## PopulationBlackberry34 -1.08199 0.24424
## PopulationBlackberry35 -0.98458 0.16551
## PopulationBlackberry36 -0.32673 0.21139
## PopulationBlackberry37 -0.61653 0.17716
## PopulationBlackberry38 -1.39482 0.33255
## PopulationBlackberry39 -0.31957 0.19152
## PopulationBlackberry40 -0.16801 0.18110
## PopulationBlackberry43 -0.57342 0.25278
## PopulationBlackberry44 0.11078 0.18110
## PopulationBlackberry45 -1.06627 0.36600
## PopulationCherry103 0.11034 0.22731
## PopulationCherry104 -1.47775 0.23736
## PopulationCherry3 -2.81107 0.28032
## PopulationCherry47 -1.04591 0.19682
## PopulationCherry50 0.71795 0.26189
## PopulationCherry51 0.88670 0.71177
## PopulationCherry52 0.39346 0.38142
## PopulationCherry6 -1.92021 0.30268
## PopulationCherry7 -0.78567 0.19891
## PopulationStrawberry42 0.08779 0.19156
## PopulationStrawberry44 -1.01054 0.20588
## PopulationStrawberry53 -0.82355 0.18506
## Test_environmentCherry 0.38655 0.21542
## Test_environmentStrawberry -0.42287 0.16978
## SA1 0.38209 0.16977
## Test_environmentBlackberry:Original_environmentCherry 0.28124 0.22182
## Test_environmentCherry:Original_environmentCherry -0.27686 0.27280
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.39178 0.29376
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## t value Pr(>|t|)
## (Intercept) 14.459 < 2e-16 ***
## PopulationBlackberry32 -13.789 < 2e-16 ***
## PopulationBlackberry33 -0.720 0.471991
## PopulationBlackberry34 -4.430 1.13e-05 ***
## PopulationBlackberry35 -5.949 4.72e-09 ***
## PopulationBlackberry36 -1.546 0.122748
## PopulationBlackberry37 -3.480 0.000539 ***
## PopulationBlackberry38 -4.194 3.17e-05 ***
## PopulationBlackberry39 -1.669 0.095747 .
## PopulationBlackberry40 -0.928 0.353952
## PopulationBlackberry43 -2.268 0.023673 *
## PopulationBlackberry44 0.612 0.540970
## PopulationBlackberry45 -2.913 0.003716 **
## PopulationCherry103 0.485 0.627572
## PopulationCherry104 -6.226 9.29e-10 ***
## PopulationCherry3 -10.028 < 2e-16 ***
## PopulationCherry47 -5.314 1.54e-07 ***
## PopulationCherry50 2.741 0.006308 **
## PopulationCherry51 1.246 0.213359
## PopulationCherry52 1.032 0.302711
## PopulationCherry6 -6.344 4.56e-10 ***
## PopulationCherry7 -3.950 8.80e-05 ***
## PopulationStrawberry42 0.458 0.646932
## PopulationStrawberry44 -4.909 1.20e-06 ***
## PopulationStrawberry53 -4.450 1.03e-05 ***
## Test_environmentCherry 1.794 0.073279 .
## Test_environmentStrawberry -2.491 0.013032 *
## SA1 2.251 0.024786 *
## Test_environmentBlackberry:Original_environmentCherry 1.268 0.205355
## Test_environmentCherry:Original_environmentCherry -1.015 0.310591
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 1.334 0.182843
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6827 on 571 degrees of freedom
## Multiple R-squared: 0.5989, Adjusted R-squared: 0.5778
## F-statistic: 28.42 on 30 and 571 DF, p-value: < 2.2e-16anova(m1)## Analysis of Variance Table
##
## Response: log(Nb_adults + 1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Population 24 323.52 13.480 28.9240 < 2.2e-16 ***
## Test_environment 2 69.66 34.828 74.7289 < 2.2e-16 ***
## SA 1 3.21 3.215 6.8980 0.008861 **
## Test_environment:Original_environment 3 0.99 0.328 0.7046 0.549525
## Residuals 571 266.12 0.466
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## F test for SA
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))## [1] 9.789562(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.05211294## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
## rsquare with two models
rsq2 <- (anova(m0)[3, 3]-anova(m1)[4, 3])/anova(m0)[3, 3]
## rsquare with single model
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))## [1] 0.6872452.3.4 Power Analysis
m2 <- lm(log(Nb_adults+1) ~ -1 + Population + Test_environment + Original_environment:Test_environment, data=data[data$Generation=="G2",])
coef(m2)
simdata <- function(seed=1, nrep = 2, npop = 3, SAincrease=0){
set.seed(seed)
## Sample populations
PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
## Create dataset
newdata <- expand.grid(Pop_new_name=1:(3*npop), Test_environment=levels(data$Test_environment), Rep=1:nrep)
newdata$Original_environment <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
## Add original names to get the right coef in the model
newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
newdata$Population <- newdata$Pop_new_name
levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
## Get design matrix
mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)
## Check that the order of the coefficients is the same
data.frame(colnames(mat), names(coef(m2))[names(coef(m2))%in%colnames(mat)])
coefsim <- coef(m2)[names(coef(m2))%in%colnames(mat)]
coefsim <- ifelse(is.na(coefsim), 0, coefsim)
## Increase fitness in sympatric environment
coefsim[grepl("PopulationBlackberry|Test_environmentCherry:Original_environmentCherry|Test_environmentStrawberry:Original_environmentStrawberry", names(coefsim))] <- coefsim[grepl("PopulationBlackberry|Test_environmentCherry:Original_environmentCherry|Test_environmentStrawberry:Original_environmentStrawberry", names(coefsim))]+SAincrease
newdata$resp <- mat%*%coefsim + rnorm(n=nrow(newdata), mean = 0, sd = sigma(m2))
m3 <- lm(resp ~ -1+Population + Test_environment + Original_environment:Test_environment, data=newdata)
data.frame(coefsim, coef(m3))
newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
m4 <- lm(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment, data=newdata)
summary(m4)
anova(m4)
## F test for SA
(Fratio_int = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
## rsquare with single model
(rsq2 <- 1-anova(m4)[4, 3]/((anova(m4)[3, 2]+anova(m4)[4, 2])/(anova(m4)[3, 1]+anova(m4)[4, 1])))
#print(rsq2)
return(pvalue_int)
}
simdata(seed=2, nrep = 30, npop = 3)
simdata(seed=2, nrep = 30, npop = 3)
simpval <- sapply(1:500, simdata, nrep = 10, npop = 3)
mean(simpval)
## Compute power
mean(simpval<0.05)
hist(simpval)
simpval10rep <- sapply(1:500, simdata, nrep = 10, npop = 3)
simpval15rep <- sapply(1:500, simdata, nrep = 15, npop = 3)
simpval20rep <- sapply(1:500, simdata, nrep = 20, npop = 3)
simpval25rep <- sapply(1:500, simdata, nrep = 25, npop = 3)
simpval30rep <- sapply(1:500, simdata, nrep = 30, npop = 3)
simpval35rep <- sapply(1:500, simdata, nrep = 35, npop = 3)
simpval40rep <- sapply(1:500, simdata, nrep = 40, npop = 3)
simpval50rep <- sapply(1:500, simdata, nrep = 50, npop = 3)
simpval60rep <- sapply(1:500, simdata, nrep = 60, npop = 3)
simpval80rep <- sapply(1:500, simdata, nrep = 80, npop = 3)
simpval100rep <- sapply(1:500, simdata, nrep = 100, npop = 3)
## Compute power
val <- data.frame(simpval10rep, simpval15rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval50rep, simpval60rep, simpval80rep, simpval100rep)
computepower <- function(x){
return(mean(x<0.05))
}
statpower <- apply(val, 2, computepower)
plot(c(10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power")2.4 Number of adults~ Number of eggs
2.4.1 Explore data
AvNb_adults <- tapply(data$Nb_adults, list(data$EggScoreFive, data$Original_environment, data$Test_environment), mean)
AvNb_adults[, , "Blackberry"][, "Cherry"]/AvNb_adults[, , "Blackberry"][, "Blackberry"]## 1 2 3 4 5
## 0.2796610 0.7280423 0.7701180 1.0788644 0.1762673AvNb_adults[, , "Cherry"][, "Blackberry"]/AvNb_adults[, , "Cherry"][, "Cherry"]## 1 2 3 4 5
## 0.0000000 0.8869732 1.1581395 0.8547641 1.5423729xyplot(Nb_adults~Nb_eggs|Original_environment*Test_environment, data=data)
m0 <- lm(log(Nb_adults+1)~log(Nb_eggs+1)+Test_environment+Original_environment, data=data)
m1 <- lm(log(Nb_adults+1)~log(Nb_eggs+1)*Test_environment+Original_environment, data=data)
m2 <- lm(log(Nb_adults+1)~log(Nb_eggs+1)*Original_environment+Test_environment, data=data)
m3 <- lm(log(Nb_adults+1)~Nb_eggs+Test_environment+Original_environment, data=data)
summary(m1)##
## Call:
## lm(formula = log(Nb_adults + 1) ~ log(Nb_eggs + 1) * Test_environment +
## Original_environment, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3193 -0.4821 0.1630 0.6116 1.7523
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -0.69352 0.50780 -1.366
## log(Nb_eggs + 1) 0.78227 0.10945 7.147
## Test_environmentCherry -0.69176 1.01112 -0.684
## Test_environmentStrawberry -0.57547 0.74678 -0.771
## Original_environmentCherry -0.29381 0.09165 -3.206
## Original_environmentStrawberry 0.15817 0.09208 1.718
## log(Nb_eggs + 1):Test_environmentCherry 0.13184 0.21267 0.620
## log(Nb_eggs + 1):Test_environmentStrawberry -0.01274 0.16191 -0.079
## Pr(>|t|)
## (Intercept) 0.17254
## log(Nb_eggs + 1) 2.61e-12 ***
## Test_environmentCherry 0.49415
## Test_environmentStrawberry 0.44125
## Original_environmentCherry 0.00142 **
## Original_environmentStrawberry 0.08635 .
## log(Nb_eggs + 1):Test_environmentCherry 0.53556
## log(Nb_eggs + 1):Test_environmentStrawberry 0.93733
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.904 on 594 degrees of freedom
## Multiple R-squared: 0.2683, Adjusted R-squared: 0.2597
## F-statistic: 31.12 on 7 and 594 DF, p-value: < 2.2e-16##Cl: density dependent larval survival rate does not seem to depend on the environment
anova(m1)## Analysis of Variance Table
##
## Response: log(Nb_adults + 1)
## Df Sum Sq Mean Sq F value Pr(>F)
## log(Nb_eggs + 1) 1 114.92 114.919 140.6112 < 2.2e-16 ***
## Test_environment 2 47.80 23.898 29.2403 7.727e-13 ***
## Original_environment 2 14.92 7.462 9.1301 0.0001243 ***
## log(Nb_eggs + 1):Test_environment 2 0.39 0.197 0.2406 0.7862571
## Residuals 594 485.47 0.817
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m2)## Analysis of Variance Table
##
## Response: log(Nb_adults + 1)
## Df Sum Sq Mean Sq F value Pr(>F)
## log(Nb_eggs + 1) 1 114.92 114.919 141.5715 < 2.2e-16 ***
## Original_environment 2 15.32 7.661 9.4372 9.231e-05 ***
## Test_environment 2 47.40 23.699 29.1953 8.050e-13 ***
## log(Nb_eggs + 1):Original_environment 2 3.69 1.843 2.2706 0.1041
## Residuals 594 482.17 0.812
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(MuMIn)
model.sel(m0, m1, m2, m3)## Model selection table
## (Int) log(Nb_egg+1) Org_env Tst_env log(Nb_egg+1):Tst_env
## m2 -1.3480 0.9262 + +
## m0 -0.7699 0.7988 + +
## m1 -0.6935 0.7823 + + +
## m3 1.9780 + +
## log(Nb_egg+1):Org_env Nb_egg df logLik AICc delta weight
## m2 + 9 -787.393 1593.1 0.00 0.521
## m0 7 -789.685 1593.6 0.47 0.412
## m1 9 -789.441 1597.2 4.10 0.067
## m3 0.008389 7 -801.139 1616.5 23.38 0.000
## Models ranked by AICc(x)2.4.2 Analysis
m1 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment+log(Nb_eggs+1), data=data[data$Generation=="G2",])
summary(m1)##
## Call:
## lm(formula = log(Nb_adults + 1) ~ Population + Test_environment +
## SA + Original_environment:Test_environment + log(Nb_eggs +
## 1), data = data[data$Generation == "G2", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.27701 -0.34196 0.05666 0.38365 1.72524
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error
## (Intercept) 0.86228 0.35909
## PopulationBlackberry32 -2.00449 0.16796
## PopulationBlackberry33 -0.14154 0.15835
## PopulationBlackberry34 -0.97556 0.23208
## PopulationBlackberry35 -0.93289 0.15714
## PopulationBlackberry36 -0.32667 0.20053
## PopulationBlackberry37 -0.57348 0.16814
## PopulationBlackberry38 -1.48328 0.31566
## PopulationBlackberry39 -0.37574 0.18182
## PopulationBlackberry40 -0.02263 0.17275
## PopulationBlackberry43 -0.40212 0.24074
## PopulationBlackberry44 0.01004 0.17225
## PopulationBlackberry45 -1.02467 0.34724
## PopulationCherry103 0.18211 0.21582
## PopulationCherry104 -1.06208 0.23104
## PopulationCherry3 -2.73259 0.26610
## PopulationCherry47 -1.09979 0.18683
## PopulationCherry50 0.61995 0.24873
## PopulationCherry51 0.86494 0.67521
## PopulationCherry52 0.44708 0.36189
## PopulationCherry6 -2.09337 0.28795
## PopulationCherry7 -0.72679 0.18883
## PopulationStrawberry42 0.02878 0.18187
## PopulationStrawberry44 -1.08884 0.19554
## PopulationStrawberry53 -0.59921 0.17776
## Test_environmentCherry 0.23651 0.20521
## Test_environmentStrawberry -0.45092 0.16110
## SA1 0.35321 0.16109
## log(Nb_eggs + 1) 0.51067 0.06358
## Test_environmentBlackberry:Original_environmentCherry 0.11866 0.21140
## Test_environmentCherry:Original_environmentCherry -0.27743 0.25879
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.26693 0.27910
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## t value Pr(>|t|)
## (Intercept) 2.401 0.016656 *
## PopulationBlackberry32 -11.934 < 2e-16 ***
## PopulationBlackberry33 -0.894 0.371793
## PopulationBlackberry34 -4.204 3.05e-05 ***
## PopulationBlackberry35 -5.937 5.06e-09 ***
## PopulationBlackberry36 -1.629 0.103865
## PopulationBlackberry37 -3.411 0.000694 ***
## PopulationBlackberry38 -4.699 3.28e-06 ***
## PopulationBlackberry39 -2.067 0.039225 *
## PopulationBlackberry40 -0.131 0.895838
## PopulationBlackberry43 -1.670 0.095398 .
## PopulationBlackberry44 0.058 0.953538
## PopulationBlackberry45 -2.951 0.003299 **
## PopulationCherry103 0.844 0.399132
## PopulationCherry104 -4.597 5.28e-06 ***
## PopulationCherry3 -10.269 < 2e-16 ***
## PopulationCherry47 -5.886 6.74e-09 ***
## PopulationCherry50 2.492 0.012970 *
## PopulationCherry51 1.281 0.200719
## PopulationCherry52 1.235 0.217181
## PopulationCherry6 -7.270 1.19e-12 ***
## PopulationCherry7 -3.849 0.000132 ***
## PopulationStrawberry42 0.158 0.874310
## PopulationStrawberry44 -5.568 3.97e-08 ***
## PopulationStrawberry53 -3.371 0.000800 ***
## Test_environmentCherry 1.153 0.249586
## Test_environmentStrawberry -2.799 0.005299 **
## SA1 2.193 0.028734 *
## log(Nb_eggs + 1) 8.032 5.54e-15 ***
## Test_environmentBlackberry:Original_environmentCherry 0.561 0.574817
## Test_environmentCherry:Original_environmentCherry -1.072 0.284151
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.956 0.339279
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6476 on 570 degrees of freedom
## Multiple R-squared: 0.6397, Adjusted R-squared: 0.6201
## F-statistic: 32.64 on 31 and 570 DF, p-value: < 2.2e-16anova(m1)## Analysis of Variance Table
##
## Response: log(Nb_adults + 1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Population 24 323.52 13.480 32.1410 < 2.2e-16 ***
## Test_environment 2 69.66 34.828 83.0404 < 2.2e-16 ***
## SA 1 3.21 3.215 7.6652 0.005813 **
## log(Nb_eggs + 1) 1 27.52 27.522 65.6206 3.338e-15 ***
## Test_environment:Original_environment 3 0.52 0.173 0.4120 0.744439
## Residuals 570 239.06 0.419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## F test for SA
(Fratio_int = (anova(m1)[3,2]/anova(m1)[5,2])/(1/anova(m1)[5,1]))## [1] 18.60447(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.02295088## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(rsq2 <- 1-anova(m1)[5, 3]/((anova(m1)[3, 2]+anova(m1)[5, 2])/(anova(m1)[3, 1]+anova(m1)[5, 1])))## [1] 0.8148531m2 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScore, data=data[data$Generation=="G2",])
summary(m2)##
## Call:
## lm(formula = log(Nb_adults + 1) ~ Population + Test_environment +
## SA + Original_environment:Test_environment + EggScore, data = data[data$Generation ==
## "G2", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3975 -0.3259 0.0477 0.3925 1.6645
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error
## (Intercept) 2.42964 0.24025
## PopulationBlackberry32 -2.07722 0.16498
## PopulationBlackberry33 -0.17843 0.15808
## PopulationBlackberry34 -1.08307 0.23151
## PopulationBlackberry35 -0.96881 0.15722
## PopulationBlackberry36 -0.37498 0.20035
## PopulationBlackberry37 -0.59598 0.16811
## PopulationBlackberry38 -1.52552 0.31662
## PopulationBlackberry39 -0.40816 0.18147
## PopulationBlackberry40 -0.08035 0.17306
## PopulationBlackberry43 -0.39924 0.24076
## PopulationBlackberry44 -0.02872 0.17257
## PopulationBlackberry45 -1.14581 0.34702
## PopulationCherry103 0.12017 0.21574
## PopulationCherry104 -1.10701 0.22934
## PopulationCherry3 -2.79795 0.26564
## PopulationCherry47 -1.14553 0.18664
## PopulationCherry50 0.56136 0.24893
## PopulationCherry51 0.81588 0.67410
## PopulationCherry52 0.37138 0.36179
## PopulationCherry6 -2.06312 0.28701
## PopulationCherry7 -0.76473 0.18843
## PopulationStrawberry42 -0.01311 0.18156
## PopulationStrawberry44 -1.10488 0.19502
## PopulationStrawberry53 -0.61495 0.17821
## Test_environmentCherry 0.19537 0.20653
## Test_environmentStrawberry -0.46864 0.16128
## SA1 0.33610 0.16217
## EggScore2 0.71482 0.12567
## EggScore3 1.00659 0.13193
## EggScore4 1.12775 0.14440
## Test_environmentBlackberry:Original_environmentCherry 0.11735 0.21231
## Test_environmentCherry:Original_environmentCherry -0.22493 0.25849
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.23700 0.28027
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## t value Pr(>|t|)
## (Intercept) 10.113 < 2e-16 ***
## PopulationBlackberry32 -12.591 < 2e-16 ***
## PopulationBlackberry33 -1.129 0.259492
## PopulationBlackberry34 -4.678 3.62e-06 ***
## PopulationBlackberry35 -6.162 1.36e-09 ***
## PopulationBlackberry36 -1.872 0.061773 .
## PopulationBlackberry37 -3.545 0.000425 ***
## PopulationBlackberry38 -4.818 1.86e-06 ***
## PopulationBlackberry39 -2.249 0.024885 *
## PopulationBlackberry40 -0.464 0.642637
## PopulationBlackberry43 -1.658 0.097822 .
## PopulationBlackberry44 -0.166 0.867882
## PopulationBlackberry45 -3.302 0.001021 **
## PopulationCherry103 0.557 0.577725
## PopulationCherry104 -4.827 1.79e-06 ***
## PopulationCherry3 -10.533 < 2e-16 ***
## PopulationCherry47 -6.138 1.57e-09 ***
## PopulationCherry50 2.255 0.024507 *
## PopulationCherry51 1.210 0.226659
## PopulationCherry52 1.027 0.305093
## PopulationCherry6 -7.188 2.08e-12 ***
## PopulationCherry7 -4.059 5.63e-05 ***
## PopulationStrawberry42 -0.072 0.942481
## PopulationStrawberry44 -5.666 2.33e-08 ***
## PopulationStrawberry53 -3.451 0.000601 ***
## Test_environmentCherry 0.946 0.344561
## Test_environmentStrawberry -2.906 0.003806 **
## SA1 2.072 0.038676 *
## EggScore2 5.688 2.06e-08 ***
## EggScore3 7.630 9.98e-14 ***
## EggScore4 7.810 2.78e-14 ***
## Test_environmentBlackberry:Original_environmentCherry 0.553 0.580680
## Test_environmentCherry:Original_environmentCherry -0.870 0.384579
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.846 0.398124
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6453 on 568 degrees of freedom
## Multiple R-squared: 0.6435, Adjusted R-squared: 0.6228
## F-statistic: 31.07 on 33 and 568 DF, p-value: < 2.2e-16anova(m2)## Analysis of Variance Table
##
## Response: log(Nb_adults + 1)
## Df Sum Sq Mean Sq F value Pr(>F)
## Population 24 323.52 13.480 32.3685 < 2.2e-16 ***
## Test_environment 2 69.66 34.828 83.6282 < 2.2e-16 ***
## SA 1 3.21 3.215 7.7195 0.005644 **
## EggScore 3 30.19 10.063 24.1633 9.951e-15 ***
## Test_environment:Original_environment 3 0.36 0.121 0.2916 0.831493
## Residuals 568 236.55 0.416
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## F test for SA
(Fratio_int = (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5,1]))## [1] 26.47414(pvalue_int = 1 - pf(Fratio_int, 1, anova(m2)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.01422731## Compute R2 = MS Interaction model without SA - MS Interaction model with SA / MS Interaction model without SA
(rsq2 <- 1-anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1])))## [1] 0.8642878## Compare with EggScoreFive
m3 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScoreFive, data=data[data$Generation=="G2",])
## Compare with EggScoreSmall
m4 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScoreSmall, data=data[data$Generation=="G2",])
## No egg information
m0 <- lm(log(Nb_adults+1) ~ Population + Test_environment + SA + Original_environment:Test_environment, data=data[data$Generation=="G2",])
model.sel(m0, m1, m2, m3, m4)## Model selection table
## (Int) Ppl SA Tst_env Org_env:Tst_env log(Nb_egg+1) EgS ESF ESS df logLik
## m4 2.3450 + + + + + 37 -567.123
## m2 2.4300 + + + + + 35 -573.037
## m3 2.4380 + + + + + 36 -572.526
## m1 0.8623 + + + + 0.5107 33 -576.218
## m0 3.2020 + + + + 32 -608.489
## AICc delta weight
## m4 1213.2 0.00 0.952
## m2 1220.5 7.29 0.025
## m3 1221.8 8.54 0.013
## m1 1222.4 9.15 0.010
## m0 1284.7 71.46 0.000
## Models ranked by AICc(x)2.4.3 Power Analysis
m2 <- lm(log(Nb_adults+1) ~ -1 + Population + Test_environment + Original_environment:Test_environment+EggScore, data=data[data$Generation=="G2",])
coef(m2)
simdata <- function(seed=1, nrep = 2, npop = 3, SAincrease=0){
set.seed(seed)
## Sample populations
PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
## Create dataset
newdata <- expand.grid(Pop_new_name=1:(3*npop), Test_environment=levels(data$Test_environment), EggScore=factor(1:4), Rep=1:nrep)
newdata$Original_environment <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
## Add original names to get the right coef in the model
newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
newdata$Population <- newdata$Pop_new_name
levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
## Get design matrix
mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment+EggScore, data=newdata)
## Check that the order of the coefficients is the same
data.frame(colnames(mat), names(coef(m2))[names(coef(m2))%in%colnames(mat)])
coefsim <- coef(m2)[names(coef(m2))%in%colnames(mat)]
coefsim <- ifelse(is.na(coefsim), 0, coefsim)
## Increase fitness in sympatric environment
coefsim[grepl("PopulationBlackberry|Test_environmentCherry:Original_environmentCherry|Test_environmentStrawberry:Original_environmentStrawberry", names(coefsim))] <- coefsim[grepl("PopulationBlackberry|Test_environmentCherry:Original_environmentCherry|Test_environmentStrawberry:Original_environmentStrawberry", names(coefsim))]+SAincrease
newdata$resp <- mat%*%coefsim + rnorm(n=nrow(newdata), mean = 0, sd = sigma(m2))
m3 <- lm(resp ~ -1+Population + Test_environment + Original_environment:Test_environment+EggScore, data=newdata)
data.frame(coefsim, coef(m3))
newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
m4 <- lm(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment+EggScore, data=newdata)
summary(m4)
anova(m4)
## F test for SA
(Fratio_int = (anova(m4)[3,2]/anova(m4)[5,2])/(1/anova(m4)[5,1]))
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)
## rsquare with single model
(rsq2 <- 1-anova(m4)[5, 3]/((anova(m4)[3, 2]+anova(m4)[5, 2])/(anova(m4)[3, 1]+anova(m4)[5, 1])))
#print(rsq2)
return(pvalue_int)
}
simdata(seed=2, nrep = 30, npop = 3)
simpval <- sapply(1:500, simdata, nrep = 10, npop = 3)
mean(simpval)
## Compute power
mean(simpval<0.05)
hist(simpval)
simpval10rep <- sapply(1:500, simdata, nrep = 10, npop = 3)
simpval15rep <- sapply(1:500, simdata, nrep = 15, npop = 3)
simpval20rep <- sapply(1:500, simdata, nrep = 20, npop = 3)
simpval25rep <- sapply(1:500, simdata, nrep = 25, npop = 3)
simpval30rep <- sapply(1:500, simdata, nrep = 30, npop = 3)
simpval35rep <- sapply(1:500, simdata, nrep = 35, npop = 3)
simpval40rep <- sapply(1:500, simdata, nrep = 40, npop = 3)
simpval50rep <- sapply(1:500, simdata, nrep = 50, npop = 3)
simpval60rep <- sapply(1:500, simdata, nrep = 60, npop = 3)
simpval80rep <- sapply(1:500, simdata, nrep = 80, npop = 3)
simpval100rep <- sapply(1:500, simdata, nrep = 100, npop = 3)
## Compute power
val <- data.frame(simpval10rep, simpval15rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval50rep, simpval60rep, simpval80rep, simpval100rep)
computepower <- function(x){
return(mean(x<0.05))
}
statpower <- apply(val, 2, computepower)
plot(c(10, 15, 20, 25, 30, 35, 40, 50, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power")2.5 Emergence rate
2.5.1 Explore data
library(lattice)
xyplot(Emergence_Rate~Nb_eggs|Original_environment*Test_environment, data=data)
xyplot(Emergence_Rate~EggScore|Original_environment*Test_environment, data=data)
bwplot(Emergence_Rate~EggScore|Original_environment*Test_environment, data=data)
bwplot(Emergence_Rate~EggScoreFive|Original_environment*Test_environment, data=data)
bwplot(Emergence_Rate~EggScoreSmall|Original_environment*Test_environment, data=data)
AvEmergenceRate <- tapply(data$Emergence_Rate, list(data$EggScoreFive, data$Original_environment, data$Test_environment), mean)
tapply(data$Emergence_Rate, list(data$EggScoreFive, data$Original_environment, data$Test_environment), length)## , , Blackberry
##
## Blackberry Cherry Strawberry
## 1 11 1 3
## 2 43 14 13
## 3 40 17 20
## 4 10 10 8
## 5 3 4 2
##
## , , Cherry
##
## Blackberry Cherry Strawberry
## 1 1 2 NA
## 2 30 10 17
## 3 50 20 21
## 4 23 11 9
## 5 2 7 NA
##
## , , Strawberry
##
## Blackberry Cherry Strawberry
## 1 7 3 6
## 2 55 22 13
## 3 38 16 16
## 4 7 6 9
## 5 NA NA 2AvEmergenceRate[, , "Blackberry"][, "Cherry"]/AvEmergenceRate[, , "Blackberry"][, "Blackberry"]## 1 2 3 4 5
## 0.2722049 0.8010658 0.7437340 0.9735309 0.1709919AvEmergenceRate[, , "Cherry"][, "Blackberry"]/AvEmergenceRate[, , "Cherry"][, "Cherry"]## 1 2 3 4 5
## 0.0000000 0.9000996 1.1452193 0.8514873 1.4918117AvEmergenceRate[, , "Blackberry"][, "Cherry"]/AvEmergenceRate[, , "Cherry"][, "Cherry"]## 1 2 3 4 5
## 0.8525097 0.9696583 1.1359488 0.9037047 0.46513152.5.2 ANR analysis
dataG2 <- data[data$Generation=="G2",]
m0 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment, data=dataG2)
dataG2$res <- residuals(m0)
m0 <- lm(res ~ -1 + Original_environment : Test_environment, data=dataG2)
summary(m0)##
## Call:
## lm(formula = res ~ -1 + Original_environment:Test_environment,
## data = dataG2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41850 -0.09344 -0.00716 0.08081 0.93906
##
## Coefficients:
## Estimate Std. Error
## Original_environmentBlackberry:Test_environmentBlackberry 0.025940 0.014200
## Original_environmentCherry:Test_environmentBlackberry -0.033385 0.021657
## Original_environmentStrawberry:Test_environmentBlackberry -0.026953 0.021657
## Original_environmentBlackberry:Test_environmentCherry -0.009440 0.014267
## Original_environmentCherry:Test_environmentCherry 0.029225 0.020773
## Original_environmentStrawberry:Test_environmentCherry -0.009799 0.021426
## Original_environmentBlackberry:Test_environmentStrawberry -0.016587 0.014200
## Original_environmentCherry:Test_environmentStrawberry 0.001584 0.021426
## Original_environmentStrawberry:Test_environmentStrawberry 0.036965 0.021657
## t value Pr(>|t|)
## Original_environmentBlackberry:Test_environmentBlackberry 1.827 0.0682 .
## Original_environmentCherry:Test_environmentBlackberry -1.541 0.1237
## Original_environmentStrawberry:Test_environmentBlackberry -1.245 0.2138
## Original_environmentBlackberry:Test_environmentCherry -0.662 0.5084
## Original_environmentCherry:Test_environmentCherry 1.407 0.1600
## Original_environmentStrawberry:Test_environmentCherry -0.457 0.6476
## Original_environmentBlackberry:Test_environmentStrawberry -1.168 0.2432
## Original_environmentCherry:Test_environmentStrawberry 0.074 0.9411
## Original_environmentStrawberry:Test_environmentStrawberry 1.707 0.0884 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1469 on 593 degrees of freedom
## Multiple R-squared: 0.02334, Adjusted R-squared: 0.008517
## F-statistic: 1.575 on 9 and 593 DF, p-value: 0.1193m0 <- lm(asin(sqrt(Emergence_Rate)) ~ -1+ Population + Test_environment + Original_environment : Test_environment, data=dataG2)
summary(m0)##
## Call:
## lm(formula = asin(sqrt(Emergence_Rate)) ~ -1 + Population + Test_environment +
## Original_environment:Test_environment, data = dataG2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41945 -0.09246 -0.00672 0.08070 0.93802
##
## Coefficients: (2 not defined because of singularities)
## Estimate Std. Error
## PopulationBlackberry31 0.63332 0.03114
## PopulationBlackberry32 0.27700 0.02673
## PopulationBlackberry33 0.57728 0.02542
## PopulationBlackberry34 0.41343 0.04653
## PopulationBlackberry35 0.43325 0.02511
## PopulationBlackberry36 0.52676 0.03841
## PopulationBlackberry37 0.52330 0.02862
## PopulationBlackberry38 0.31701 0.06781
## PopulationBlackberry39 0.51650 0.03277
## PopulationBlackberry40 0.66341 0.02978
## PopulationBlackberry43 0.63278 0.04895
## PopulationBlackberry44 0.60013 0.02978
## PopulationBlackberry45 0.40441 0.07536
## PopulationCherry103 0.69122 0.04398
## PopulationCherry104 0.42678 0.04646
## PopulationCherry3 0.16504 0.05681
## PopulationCherry47 0.38239 0.03622
## PopulationCherry50 0.77006 0.05241
## PopulationCherry51 0.92529 0.15428
## PopulationCherry52 0.78469 0.08027
## PopulationCherry6 0.24721 0.06208
## PopulationCherry7 0.45095 0.03677
## PopulationStrawberry42 0.67561 0.03518
## PopulationStrawberry44 0.42915 0.03860
## PopulationStrawberry53 0.57813 0.03315
## Test_environmentCherry -0.07930 0.02054
## Test_environmentStrawberry -0.18264 0.02049
## Test_environmentBlackberry:Original_environmentCherry -0.07859 0.03741
## Test_environmentCherry:Original_environmentCherry 0.02052 0.03673
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry -0.10655 0.03733
## Test_environmentCherry:Original_environmentStrawberry -0.05393 0.03722
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## t value Pr(>|t|)
## PopulationBlackberry31 20.335 < 2e-16 ***
## PopulationBlackberry32 10.362 < 2e-16 ***
## PopulationBlackberry33 22.709 < 2e-16 ***
## PopulationBlackberry34 8.885 < 2e-16 ***
## PopulationBlackberry35 17.257 < 2e-16 ***
## PopulationBlackberry36 13.715 < 2e-16 ***
## PopulationBlackberry37 18.285 < 2e-16 ***
## PopulationBlackberry38 4.675 3.68e-06 ***
## PopulationBlackberry39 15.763 < 2e-16 ***
## PopulationBlackberry40 22.276 < 2e-16 ***
## PopulationBlackberry43 12.927 < 2e-16 ***
## PopulationBlackberry44 20.151 < 2e-16 ***
## PopulationBlackberry45 5.366 1.17e-07 ***
## PopulationCherry103 15.718 < 2e-16 ***
## PopulationCherry104 9.187 < 2e-16 ***
## PopulationCherry3 2.905 0.003814 **
## PopulationCherry47 10.556 < 2e-16 ***
## PopulationCherry50 14.693 < 2e-16 ***
## PopulationCherry51 5.997 3.56e-09 ***
## PopulationCherry52 9.776 < 2e-16 ***
## PopulationCherry6 3.982 7.72e-05 ***
## PopulationCherry7 12.265 < 2e-16 ***
## PopulationStrawberry42 19.205 < 2e-16 ***
## PopulationStrawberry44 11.117 < 2e-16 ***
## PopulationStrawberry53 17.439 < 2e-16 ***
## Test_environmentCherry -3.861 0.000126 ***
## Test_environmentStrawberry -8.914 < 2e-16 ***
## Test_environmentBlackberry:Original_environmentCherry -2.101 0.036095 *
## Test_environmentCherry:Original_environmentCherry 0.559 0.576564
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry -2.854 0.004474 **
## Test_environmentCherry:Original_environmentStrawberry -1.449 0.147919
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1497 on 571 degrees of freedom
## Multiple R-squared: 0.9001, Adjusted R-squared: 0.8947
## F-statistic: 166 on 31 and 571 DF, p-value: < 2.2e-16emmeans::emmeans(m0, list(pairwise ~ Original_environment:Test_environment), adjust = "tukey") ### Warning in (function (object, at, cov.reduce = mean, cov.keep = get_emm_option("cov.keep"), : There are unevaluated constants in the response formula
## Auto-detection of the response transformation may be incorrect## NOTE: A nesting structure was detected in the fitted model:
## Population %in% (Test_environment*Original_environment), Original_environment %in% Test_environment## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $`emmeans of Original_environment, Test_environment`
## Original_environment Test_environment emmean SE df lower.CL upper.CL
## Blackberry Blackberry 0.501 0.0161 571 0.470 0.533
## Cherry Blackberry 0.460 0.0273 571 0.406 0.513
## Strawberry Blackberry 0.454 0.0224 571 0.410 0.498
## Blackberry Cherry 0.422 0.0163 571 0.390 0.454
## Cherry Cherry 0.431 0.0239 571 0.384 0.478
## Strawberry Cherry 0.428 0.0222 571 0.384 0.471
## Blackberry Strawberry 0.319 0.0165 571 0.286 0.351
## Cherry Strawberry 0.307 0.0242 571 0.260 0.355
## Strawberry Strawberry 0.378 0.0224 571 0.334 0.422
##
## Results are averaged over the levels of: Population
## Results are given on the asin(sqrt(mu)) (not the response) scale.
## Confidence level used: 0.95
##
## $`pairwise differences of Original_environment, Test_environment`
## contrast estimate SE df t.ratio
## Blackberry Blackberry - Cherry Blackberry 0.04184 0.0317 571 1.319
## Blackberry Blackberry - Strawberry Blackberry 0.04702 0.0276 571 1.704
## Blackberry Blackberry - Blackberry Cherry 0.07930 0.0205 571 3.861
## Blackberry Blackberry - Cherry Cherry 0.07042 0.0289 571 2.441
## Blackberry Blackberry - Strawberry Cherry 0.07370 0.0274 571 2.688
## Blackberry Blackberry - Blackberry Strawberry 0.18264 0.0205 571 8.914
## Blackberry Blackberry - Cherry Strawberry 0.19428 0.0291 571 6.683
## Blackberry Blackberry - Strawberry Strawberry 0.12311 0.0276 571 4.460
## Cherry Blackberry - Strawberry Blackberry 0.00518 0.0353 571 0.147
## Cherry Blackberry - Blackberry Cherry 0.03746 0.0318 571 1.177
## Cherry Blackberry - Cherry Cherry 0.02858 0.0336 571 0.852
## Cherry Blackberry - Strawberry Cherry 0.03185 0.0352 571 0.905
## Cherry Blackberry - Blackberry Strawberry 0.14080 0.0319 571 4.414
## Cherry Blackberry - Cherry Strawberry 0.15244 0.0340 571 4.488
## Cherry Blackberry - Strawberry Strawberry 0.08127 0.0353 571 2.300
## Strawberry Blackberry - Blackberry Cherry 0.03228 0.0277 571 1.164
## Strawberry Blackberry - Cherry Cherry 0.02340 0.0328 571 0.713
## Strawberry Blackberry - Strawberry Cherry 0.02667 0.0310 571 0.859
## Strawberry Blackberry - Blackberry Strawberry 0.13562 0.0278 571 4.875
## Strawberry Blackberry - Cherry Strawberry 0.14726 0.0330 571 4.464
## Strawberry Blackberry - Strawberry Strawberry 0.07609 0.0312 571 2.438
## Blackberry Cherry - Cherry Cherry -0.00888 0.0290 571 -0.306
## Blackberry Cherry - Strawberry Cherry -0.00560 0.0275 571 -0.203
## Blackberry Cherry - Blackberry Strawberry 0.10334 0.0205 571 5.033
## Blackberry Cherry - Cherry Strawberry 0.11498 0.0292 571 3.938
## Blackberry Cherry - Strawberry Strawberry 0.04381 0.0277 571 1.580
## Cherry Cherry - Strawberry Cherry 0.00328 0.0326 571 0.100
## Cherry Cherry - Blackberry Strawberry 0.11222 0.0291 571 3.862
## Cherry Cherry - Cherry Strawberry 0.12386 0.0304 571 4.068
## Cherry Cherry - Strawberry Strawberry 0.05269 0.0328 571 1.607
## Strawberry Cherry - Blackberry Strawberry 0.10895 0.0276 571 3.943
## Strawberry Cherry - Cherry Strawberry 0.12059 0.0328 571 3.673
## Strawberry Cherry - Strawberry Strawberry 0.04941 0.0310 571 1.592
## Blackberry Strawberry - Cherry Strawberry 0.01164 0.0293 571 0.398
## Blackberry Strawberry - Strawberry Strawberry -0.05953 0.0278 571 -2.140
## Cherry Strawberry - Strawberry Strawberry -0.07117 0.0330 571 -2.158
## p.value
## 0.9254
## 0.7443
## 0.0040
## 0.2640
## 0.1543
## <.0001
## <.0001
## 0.0003
## 1.0000
## 0.9611
## 0.9952
## 0.9927
## 0.0004
## 0.0003
## 0.3440
## 0.9636
## 0.9986
## 0.9948
## <.0001
## 0.0003
## 0.2655
## 1.0000
## 1.0000
## <.0001
## 0.0030
## 0.8157
## 1.0000
## 0.0040
## 0.0018
## 0.8011
## 0.0029
## 0.0080
## 0.8093
## 1.0000
## 0.4469
## 0.4352
##
## Results are averaged over the levels of: Population
## Note: contrasts are still on the asin.sqrt scale
## P value adjustment: tukey method for comparing a family of 9 estimatesdataG2$Nb_adults_corr <- dataG2$Nb_adults
dataG2$Nb_adults_corr[dataG2$Nb_adults>dataG2$Nb_eggs] <- dataG2$Nb_eggs[dataG2$Nb_adults>dataG2$Nb_eggs]
m0 <- glm(cbind(Nb_adults_corr, Nb_eggs-Nb_adults_corr) ~ -1+ Population + Test_environment, family=binomial, data=dataG2)
hist(residuals(m0))
m0 <- lm(residuals(m0) ~ -1 + Original_environment : Test_environment, data=dataG2)
summary(m0)##
## Call:
## lm(formula = residuals(m0) ~ -1 + Original_environment:Test_environment,
## data = dataG2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.3245 -1.8429 -0.1601 1.6310 12.7946
##
## Coefficients:
## Estimate Std. Error
## Original_environmentBlackberry:Test_environmentBlackberry 0.2566 0.2570
## Original_environmentCherry:Test_environmentBlackberry -0.7885 0.3920
## Original_environmentStrawberry:Test_environmentBlackberry -0.4115 0.3920
## Original_environmentBlackberry:Test_environmentCherry -0.2386 0.2582
## Original_environmentCherry:Test_environmentCherry 0.5971 0.3760
## Original_environmentStrawberry:Test_environmentCherry 0.0459 0.3878
## Original_environmentBlackberry:Test_environmentStrawberry -0.5650 0.2570
## Original_environmentCherry:Test_environmentStrawberry -0.2881 0.3878
## Original_environmentStrawberry:Test_environmentStrawberry 0.4425 0.3920
## t value Pr(>|t|)
## Original_environmentBlackberry:Test_environmentBlackberry 0.998 0.3186
## Original_environmentCherry:Test_environmentBlackberry -2.011 0.0447 *
## Original_environmentStrawberry:Test_environmentBlackberry -1.050 0.2943
## Original_environmentBlackberry:Test_environmentCherry -0.924 0.3559
## Original_environmentCherry:Test_environmentCherry 1.588 0.1129
## Original_environmentStrawberry:Test_environmentCherry 0.118 0.9058
## Original_environmentBlackberry:Test_environmentStrawberry -2.198 0.0283 *
## Original_environmentCherry:Test_environmentStrawberry -0.743 0.4579
## Original_environmentStrawberry:Test_environmentStrawberry 1.129 0.2594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.659 on 593 degrees of freedom
## Multiple R-squared: 0.02658, Adjusted R-squared: 0.0118
## F-statistic: 1.799 on 9 and 593 DF, p-value: 0.06551emmeans::emmeans(m0, list(pairwise ~ Original_environment:Test_environment), adjust = "tukey") ### NOTE: A nesting structure was detected in the fitted model:
## Original_environment %in% (Original_environment*Test_environment), Test_environment %in% (Original_environment*Test_environment)
## Note: Use 'contrast(regrid(object), ...)' to obtain contrasts of back-transformed estimates## $`emmeans of Original_environment, Test_environment`
## Original_environment Test_environment emmean SE df lower.CL upper.CL
## Blackberry Blackberry 0.2566 0.257 593 -0.248 0.7614
## Cherry Blackberry -0.7885 0.392 593 -1.558 -0.0185
## Strawberry Blackberry -0.4115 0.392 593 -1.181 0.3585
## Blackberry Cherry -0.2386 0.258 593 -0.746 0.2686
## Cherry Cherry 0.5971 0.376 593 -0.141 1.3355
## Strawberry Cherry 0.0459 0.388 593 -0.716 0.8076
## Blackberry Strawberry -0.5650 0.257 593 -1.070 -0.0602
## Cherry Strawberry -0.2881 0.388 593 -1.050 0.4736
## Strawberry Strawberry 0.4425 0.392 593 -0.327 1.2125
##
## Results are given on the residuals (not the response) scale.
## Confidence level used: 0.95
##
## $`pairwise differences of Original_environment, Test_environment`
## contrast estimate SE df t.ratio
## Blackberry Blackberry - Cherry Blackberry 1.0450 0.469 593 2.229
## Blackberry Blackberry - Strawberry Blackberry 0.6680 0.469 593 1.425
## Blackberry Blackberry - Blackberry Cherry 0.4952 0.364 593 1.359
## Blackberry Blackberry - Cherry Cherry -0.3405 0.455 593 -0.748
## Blackberry Blackberry - Strawberry Cherry 0.2107 0.465 593 0.453
## Blackberry Blackberry - Blackberry Strawberry 0.8216 0.364 593 2.260
## Blackberry Blackberry - Cherry Strawberry 0.5447 0.465 593 1.171
## Blackberry Blackberry - Strawberry Strawberry -0.1860 0.469 593 -0.397
## Cherry Blackberry - Strawberry Blackberry -0.3770 0.554 593 -0.680
## Cherry Blackberry - Blackberry Cherry -0.5499 0.469 593 -1.171
## Cherry Blackberry - Cherry Cherry -1.3855 0.543 593 -2.551
## Cherry Blackberry - Strawberry Cherry -0.8344 0.551 593 -1.513
## Cherry Blackberry - Blackberry Strawberry -0.2235 0.469 593 -0.477
## Cherry Blackberry - Cherry Strawberry -0.5004 0.551 593 -0.907
## Cherry Blackberry - Strawberry Strawberry -1.2310 0.554 593 -2.220
## Strawberry Blackberry - Blackberry Cherry -0.1729 0.469 593 -0.368
## Strawberry Blackberry - Cherry Cherry -1.0085 0.543 593 -1.857
## Strawberry Blackberry - Strawberry Cherry -0.4574 0.551 593 -0.829
## Strawberry Blackberry - Blackberry Strawberry 0.1535 0.469 593 0.327
## Strawberry Blackberry - Cherry Strawberry -0.1234 0.551 593 -0.224
## Strawberry Blackberry - Strawberry Strawberry -0.8540 0.554 593 -1.540
## Blackberry Cherry - Cherry Cherry -0.8356 0.456 593 -1.832
## Blackberry Cherry - Strawberry Cherry -0.2845 0.466 593 -0.611
## Blackberry Cherry - Blackberry Strawberry 0.3264 0.364 593 0.896
## Blackberry Cherry - Cherry Strawberry 0.0495 0.466 593 0.106
## Blackberry Cherry - Strawberry Strawberry -0.6811 0.469 593 -1.451
## Cherry Cherry - Strawberry Cherry 0.5512 0.540 593 1.020
## Cherry Cherry - Blackberry Strawberry 1.1620 0.455 593 2.551
## Cherry Cherry - Cherry Strawberry 0.8851 0.540 593 1.639
## Cherry Cherry - Strawberry Strawberry 0.1545 0.543 593 0.284
## Strawberry Cherry - Blackberry Strawberry 0.6109 0.465 593 1.313
## Strawberry Cherry - Cherry Strawberry 0.3340 0.548 593 0.609
## Strawberry Cherry - Strawberry Strawberry -0.3966 0.551 593 -0.719
## Blackberry Strawberry - Cherry Strawberry -0.2769 0.465 593 -0.595
## Blackberry Strawberry - Strawberry Strawberry -1.0075 0.469 593 -2.149
## Cherry Strawberry - Strawberry Strawberry -0.7306 0.551 593 -1.325
## p.value
## 0.3879
## 0.8880
## 0.9125
## 0.9981
## 1.0000
## 0.3684
## 0.9623
## 1.0000
## 0.9990
## 0.9622
## 0.2102
## 0.8493
## 0.9999
## 0.9925
## 0.3937
## 1.0000
## 0.6440
## 0.9960
## 1.0000
## 1.0000
## 0.8360
## 0.6608
## 0.9996
## 0.9932
## 1.0000
## 0.8773
## 0.9839
## 0.2099
## 0.7831
## 1.0000
## 0.9274
## 0.9996
## 0.9985
## 0.9996
## 0.4407
## 0.9237
##
## Note: contrasts are still on the residuals scale
## P value adjustment: tukey method for comparing a family of 9 estimates2.5.3 Explore data with blue
m0 <- lm(Emergence_Rate~Original_environment*Test_environment, data=data[data$Generation=="G2",])
summary(m0)##
## Call:
## lm(formula = Emergence_Rate ~ Original_environment * Test_environment,
## data = data[data$Generation == "G2", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.26484 -0.09978 -0.02467 0.07237 0.73516
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 0.26484 0.01330
## Original_environmentCherry -0.07926 0.02426
## Original_environmentStrawberry -0.04419 0.02426
## Test_environmentCherry -0.07253 0.01886
## Test_environmentStrawberry -0.13033 0.01881
## Original_environmentCherry:Test_environmentCherry 0.06711 0.03385
## Original_environmentStrawberry:Test_environmentCherry 0.04696 0.03421
## Original_environmentCherry:Test_environmentStrawberry 0.05509 0.03418
## Original_environmentStrawberry:Test_environmentStrawberry 0.07162 0.03431
## t value Pr(>|t|)
## (Intercept) 19.909 < 2e-16 ***
## Original_environmentCherry -3.267 0.001150 **
## Original_environmentStrawberry -1.822 0.069020 .
## Test_environmentCherry -3.846 0.000133 ***
## Test_environmentStrawberry -6.928 1.12e-11 ***
## Original_environmentCherry:Test_environmentCherry 1.982 0.047891 *
## Original_environmentStrawberry:Test_environmentCherry 1.373 0.170346
## Original_environmentCherry:Test_environmentStrawberry 1.612 0.107534
## Original_environmentStrawberry:Test_environmentStrawberry 2.088 0.037265 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1376 on 593 degrees of freedom
## Multiple R-squared: 0.1044, Adjusted R-squared: 0.09231
## F-statistic: 8.64 on 8 and 593 DF, p-value: 3.499e-11tapply(data$Emergence_Rate[data$Generation=="G2"], list(data$Original_environment[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)## Blackberry Cherry Strawberry
## Blackberry 0.2648444 0.1923140 0.1345134
## Cherry 0.1855892 0.1801664 0.1103531
## Strawberry 0.2206513 0.1950778 0.1619433tapply(data$Emergence_Rate[data$Generation=="G2"], list(data$Original_environment[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), var)## Blackberry Cherry Strawberry
## Blackberry 0.03491514 0.013749525 0.01436498
## Cherry 0.02922031 0.019467766 0.01399160
## Strawberry 0.01555929 0.009464032 0.01139614range(data$Nb_adults)## [1] 0 108## Check number of eggs and adults
tapply(data$Nb_eggs[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)## Blackberry Cherry Strawberry
## Blackberry31 106.44444 130.33333 118.88889
## Blackberry32 61.23077 84.46154 74.00000
## Blackberry33 125.13333 141.28571 95.66667
## Blackberry34 105.75000 95.66667 83.50000
## Blackberry35 93.60000 123.73333 96.31250
## Blackberry36 118.00000 144.00000 92.83333
## Blackberry37 108.09091 114.45455 101.27273
## Blackberry38 125.50000 141.00000 143.00000
## Blackberry39 133.75000 146.00000 117.12500
## Blackberry40 88.30000 92.10000 80.40000
## Blackberry43 69.33333 90.00000 86.33333
## Blackberry44 129.20000 155.70000 129.70000
## Blackberry45 92.00000 113.00000 110.00000
## Cherry103 129.66667 118.00000 89.83333
## Cherry104 59.80000 61.83333 55.20000
## Cherry3 109.66667 143.66667 78.00000
## Cherry47 142.00000 162.71429 128.23077
## Cherry50 148.33333 185.00000 126.00000
## Cherry51 139.00000 NA NA
## Cherry52 129.50000 96.00000 103.00000
## Cherry6 226.50000 203.50000 124.00000
## Cherry7 137.75000 118.69231 97.50000
## Strawberry42 140.46667 131.81250 146.66667
## Strawberry44 150.60000 142.30000 139.70000
## Strawberry53 89.00000 99.23810 73.09524tapply(data$Nb_adults[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean)## Blackberry Cherry Strawberry
## Blackberry31 38.444444 36.777778 19.2222222
## Blackberry32 1.923077 6.307692 0.6923077
## Blackberry33 30.800000 34.000000 16.6666667
## Blackberry34 13.500000 17.666667 8.0000000
## Blackberry35 14.533333 23.533333 5.2500000
## Blackberry36 28.800000 23.166667 14.0000000
## Blackberry37 34.727273 24.090909 8.2727273
## Blackberry38 15.500000 7.000000 1.0000000
## Blackberry39 28.375000 26.750000 19.0000000
## Blackberry40 38.700000 17.300000 23.9000000
## Blackberry43 45.000000 13.750000 11.6666667
## Blackberry44 48.900000 27.300000 25.6000000
## Blackberry45 11.500000 30.000000 2.0000000
## Cherry103 35.666667 41.857143 20.8333333
## Cherry104 8.200000 10.833333 3.0000000
## Cherry3 1.000000 3.000000 0.0000000
## Cherry47 15.833333 18.571429 5.9230769
## Cherry50 63.000000 67.750000 46.0000000
## Cherry51 78.000000 NA NA
## Cherry52 65.500000 39.000000 18.0000000
## Cherry6 5.500000 10.000000 1.3333333
## Cherry7 16.416667 16.076923 9.9166667
## Strawberry42 45.733333 35.000000 26.3333333
## Strawberry44 13.100000 13.400000 9.9000000
## Strawberry53 17.285714 17.000000 11.6190476## Check for the presence of negative correlations
m1 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment, data=data[data$Generation=="G2",])
data$resid <- residuals(m1)
meanbypopbytestenv <- as.data.frame(tapply(data$resid[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean))
## Cherry ~ Blackberry
plot(meanbypopbytestenv$Blackberry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Cherry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Strawberry ~ Cherry
plot(meanbypopbytestenv$Cherry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Cherry", ylab="Strawberry", pch=16)
legend("bottomright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Strawberry ~ Blackberry
plot(meanbypopbytestenv$Blackberry~ meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Strawberry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
2.5.4 Explore data with BLUPS
m0 <- lmer(asin(sqrt(Emergence_Rate))~ Population + Test_environment + (0 + Test_environment | Population), data = data[data$Generation=="G2",])## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## unable to evaluate scaled gradient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge: degenerate Hessian with 2 negative eigenvalues## Correlation between fixed effects
plot(coef(m1), fixef(m0))
summary(m0)## Linear mixed model fit by REML ['lmerMod']
## Formula: asin(sqrt(Emergence_Rate)) ~ Population + Test_environment +
## (0 + Test_environment | Population)
## Data: data[data$Generation == "G2", ]
##
## REML criterion at convergence: -480.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4905 -0.6208 -0.0492 0.5313 5.5124
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## Population Test_environmentBlackberry 0.0009946 0.03154
## Test_environmentCherry 0.0114799 0.10714 0.08
## Test_environmentStrawberry 0.0080802 0.08989 0.33 0.79
## Residual 0.0206939 0.14385
## Number of obs: 602, groups: Population, 25
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.65326 0.05211 12.536
## PopulationBlackberry32 -0.42343 0.07005 -6.045
## PopulationBlackberry33 -0.11104 0.06898 -1.610
## PopulationBlackberry34 -0.31902 0.08364 -3.814
## PopulationBlackberry35 -0.25882 0.06897 -3.753
## PopulationBlackberry36 -0.14211 0.07997 -1.777
## PopulationBlackberry37 -0.08778 0.07144 -1.229
## PopulationBlackberry38 -0.33401 0.09810 -3.405
## PopulationBlackberry39 -0.17511 0.07451 -2.350
## PopulationBlackberry40 0.03641 0.07230 0.504
## PopulationBlackberry43 0.11160 0.08780 1.271
## PopulationBlackberry44 -0.02326 0.07230 -0.322
## PopulationBlackberry45 -0.27941 0.10044 -2.782
## PopulationCherry103 -0.02717 0.07774 -0.350
## PopulationCherry104 -0.27661 0.08008 -3.454
## PopulationCherry3 -0.55496 0.08815 -6.296
## PopulationCherry47 -0.32855 0.07066 -4.649
## PopulationCherry50 0.05616 0.08738 0.643
## PopulationCherry51 0.19344 0.15622 1.238
## PopulationCherry52 0.09981 0.10044 0.994
## PopulationCherry6 -0.47735 0.09491 -5.030
## PopulationCherry7 -0.28451 0.07068 -4.025
## PopulationStrawberry42 -0.04642 0.06896 -0.673
## PopulationStrawberry44 -0.34414 0.07230 -4.760
## PopulationStrawberry53 -0.19139 0.06681 -2.865
## Test_environmentCherry -0.04937 0.02799 -1.764
## Test_environmentStrawberry -0.14711 0.02380 -6.180##
## Correlation matrix not shown by default, as p = 27 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 2 negative eigenvaluesBLUP <- ranef(m0)$Population
pairs(BLUP)
## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentBlackberry, meanbypopbytestenv$Blackberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentCherry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentStrawberry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("bottomright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
2.5.5 Analysis
data$SA <- as.factor(ifelse(data$Original_environment==data$Test_environment, 1, 0))
m1 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment, data=data[data$Generation=="G2",])
summary(m1)##
## Call:
## lm(formula = asin(sqrt(Emergence_Rate)) ~ Population + Test_environment +
## SA + Original_environment:Test_environment, data = data[data$Generation ==
## "G2", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41945 -0.09246 -0.00672 0.08070 0.93802
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error
## (Intercept) 0.5793941 0.0485461
## PopulationBlackberry32 -0.3563256 0.0374731
## PopulationBlackberry33 -0.0560407 0.0365927
## PopulationBlackberry34 -0.2198934 0.0535505
## PopulationBlackberry35 -0.2000730 0.0362889
## PopulationBlackberry36 -0.1065683 0.0463480
## PopulationBlackberry37 -0.1100221 0.0388418
## PopulationBlackberry38 -0.3163117 0.0729120
## PopulationBlackberry39 -0.1168236 0.0419913
## PopulationBlackberry40 0.0300903 0.0397061
## PopulationBlackberry43 -0.0005464 0.0554218
## PopulationBlackberry44 -0.0331930 0.0397061
## PopulationBlackberry45 -0.2289148 0.0802466
## PopulationCherry103 0.0578950 0.0498387
## PopulationCherry104 -0.2065488 0.0520413
## PopulationCherry3 -0.4682877 0.0614604
## PopulationCherry47 -0.2509357 0.0431538
## PopulationCherry50 0.1367322 0.0574191
## PopulationCherry51 0.2919698 0.1560561
## PopulationCherry52 0.1513627 0.0836261
## PopulationCherry6 -0.3861147 0.0663642
## PopulationCherry7 -0.1823742 0.0436110
## PopulationStrawberry42 -0.0116455 0.0419997
## PopulationStrawberry44 -0.2581011 0.0451386
## PopulationStrawberry53 -0.1091287 0.0405748
## Test_environmentCherry -0.0253687 0.0472313
## Test_environmentStrawberry -0.1287121 0.0372242
## SA1 0.0539301 0.0372219
## Test_environmentBlackberry:Original_environmentCherry -0.0246612 0.0486347
## Test_environmentCherry:Original_environmentCherry -0.0334106 0.0598113
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.0013057 0.0644075
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## t value Pr(>|t|)
## (Intercept) 11.935 < 2e-16 ***
## PopulationBlackberry32 -9.509 < 2e-16 ***
## PopulationBlackberry33 -1.531 0.126207
## PopulationBlackberry34 -4.106 4.61e-05 ***
## PopulationBlackberry35 -5.513 5.34e-08 ***
## PopulationBlackberry36 -2.299 0.021847 *
## PopulationBlackberry37 -2.833 0.004780 **
## PopulationBlackberry38 -4.338 1.70e-05 ***
## PopulationBlackberry39 -2.782 0.005580 **
## PopulationBlackberry40 0.758 0.448868
## PopulationBlackberry43 -0.010 0.992137
## PopulationBlackberry44 -0.836 0.403523
## PopulationBlackberry45 -2.853 0.004493 **
## PopulationCherry103 1.162 0.245864
## PopulationCherry104 -3.969 8.14e-05 ***
## PopulationCherry3 -7.619 1.07e-13 ***
## PopulationCherry47 -5.815 1.01e-08 ***
## PopulationCherry50 2.381 0.017578 *
## PopulationCherry51 1.871 0.061866 .
## PopulationCherry52 1.810 0.070822 .
## PopulationCherry6 -5.818 9.92e-09 ***
## PopulationCherry7 -4.182 3.35e-05 ***
## PopulationStrawberry42 -0.277 0.781668
## PopulationStrawberry44 -5.718 1.74e-08 ***
## PopulationStrawberry53 -2.690 0.007364 **
## Test_environmentCherry -0.537 0.591397
## Test_environmentStrawberry -3.458 0.000585 ***
## SA1 1.449 0.147919
## Test_environmentBlackberry:Original_environmentCherry -0.507 0.612301
## Test_environmentCherry:Original_environmentCherry -0.559 0.576653
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.020 0.983833
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1497 on 571 degrees of freedom
## Multiple R-squared: 0.493, Adjusted R-squared: 0.4664
## F-statistic: 18.51 on 30 and 571 DF, p-value: < 2.2e-16anova(m1)## Analysis of Variance Table
##
## Response: asin(sqrt(Emergence_Rate))
## Df Sum Sq Mean Sq F value Pr(>F)
## Population 24 10.0931 0.42054 18.7710 < 2.2e-16 ***
## Test_environment 2 2.0401 1.02004 45.5294 < 2.2e-16 ***
## SA 1 0.2890 0.28899 12.8989 0.0003571 ***
## Test_environment:Original_environment 3 0.0186 0.00621 0.2773 0.8417788
## Residuals 571 12.7926 0.02240
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1## F test for SA
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))## [1] 46.51257(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.006448859## Compute Rsquare
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))## [1] 0.9192124## Correct for number of eggs
m2 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment+log(Nb_eggs+1), data=data[data$Generation=="G2",])
## F test for SA
(Fratio_int = (anova(m2)[3,2]/anova(m2)[5,2])/(1/anova(m2)[5,1]))## [1] 108.7824(pvalue_int = 1 - pf(Fratio_int, 1, anova(m2)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.001881239## Compute Rsquare
(rsq2 <- 1-anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1])))## [1] 0.9642162m3 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScore, data=data[data$Generation=="G2",])
summary(m3)##
## Call:
## lm(formula = asin(sqrt(Emergence_Rate)) ~ Population + Test_environment +
## SA + Original_environment:Test_environment + EggScore, data = data[data$Generation ==
## "G2", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.43011 -0.09039 -0.00357 0.07908 0.93919
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error
## (Intercept) 0.602001 0.055208
## PopulationBlackberry32 -0.377760 0.037911
## PopulationBlackberry33 -0.055998 0.036326
## PopulationBlackberry34 -0.229305 0.053200
## PopulationBlackberry35 -0.210254 0.036129
## PopulationBlackberry36 -0.107246 0.046039
## PopulationBlackberry37 -0.120148 0.038631
## PopulationBlackberry38 -0.321473 0.072756
## PopulationBlackberry39 -0.108157 0.041700
## PopulationBlackberry40 0.011493 0.039768
## PopulationBlackberry43 -0.023178 0.055325
## PopulationBlackberry44 -0.027083 0.039655
## PopulationBlackberry45 -0.237178 0.079741
## PopulationCherry103 0.045284 0.049575
## PopulationCherry104 -0.240476 0.052701
## PopulationCherry3 -0.480668 0.061042
## PopulationCherry47 -0.240715 0.042888
## PopulationCherry50 0.158899 0.057202
## PopulationCherry51 0.278850 0.154902
## PopulationCherry52 0.136201 0.083137
## PopulationCherry6 -0.367035 0.065953
## PopulationCherry7 -0.191145 0.043298
## PopulationStrawberry42 -0.008004 0.041722
## PopulationStrawberry44 -0.255111 0.044814
## PopulationStrawberry53 -0.136834 0.040950
## Test_environmentCherry 0.001684 0.047459
## Test_environmentStrawberry -0.117301 0.037060
## SA1 0.069153 0.037266
## EggScore2 -0.016370 0.028877
## EggScore3 -0.036638 0.030316
## EggScore4 -0.088070 0.033183
## Test_environmentBlackberry:Original_environmentCherry 0.002489 0.048787
## Test_environmentCherry:Original_environmentCherry -0.044893 0.059398
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.029211 0.064404
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## t value Pr(>|t|)
## (Intercept) 10.904 < 2e-16 ***
## PopulationBlackberry32 -9.964 < 2e-16 ***
## PopulationBlackberry33 -1.542 0.123741
## PopulationBlackberry34 -4.310 1.92e-05 ***
## PopulationBlackberry35 -5.820 9.86e-09 ***
## PopulationBlackberry36 -2.329 0.020184 *
## PopulationBlackberry37 -3.110 0.001964 **
## PopulationBlackberry38 -4.419 1.19e-05 ***
## PopulationBlackberry39 -2.594 0.009741 **
## PopulationBlackberry40 0.289 0.772692
## PopulationBlackberry43 -0.419 0.675418
## PopulationBlackberry44 -0.683 0.494907
## PopulationBlackberry45 -2.974 0.003061 **
## PopulationCherry103 0.913 0.361395
## PopulationCherry104 -4.563 6.18e-06 ***
## PopulationCherry3 -7.874 1.75e-14 ***
## PopulationCherry47 -5.613 3.12e-08 ***
## PopulationCherry50 2.778 0.005653 **
## PopulationCherry51 1.800 0.072364 .
## PopulationCherry52 1.638 0.101919
## PopulationCherry6 -5.565 4.04e-08 ***
## PopulationCherry7 -4.415 1.21e-05 ***
## PopulationStrawberry42 -0.192 0.847934
## PopulationStrawberry44 -5.693 2.01e-08 ***
## PopulationStrawberry53 -3.341 0.000888 ***
## Test_environmentCherry 0.035 0.971714
## Test_environmentStrawberry -3.165 0.001633 **
## SA1 1.856 0.064023 .
## EggScore2 -0.567 0.571005
## EggScore3 -1.209 0.227343
## EggScore4 -2.654 0.008175 **
## Test_environmentBlackberry:Original_environmentCherry 0.051 0.959323
## Test_environmentCherry:Original_environmentCherry -0.756 0.450084
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.454 0.650315
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1483 on 568 degrees of freedom
## Multiple R-squared: 0.505, Adjusted R-squared: 0.4762
## F-statistic: 17.56 on 33 and 568 DF, p-value: < 2.2e-16## F test for SA
(Fratio_int = (anova(m3)[3,2]/anova(m3)[5,2])/(1/anova(m3)[5,1]))## [1] 64.43531(pvalue_int = 1 - pf(Fratio_int, 1, anova(m3)[5,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.004036806## Compute Rsquare
(rsq2 <- 1-anova(m2)[5, 3]/((anova(m2)[3, 2]+anova(m2)[5, 2])/(anova(m2)[3, 1]+anova(m2)[5, 1])))## [1] 0.9642162## Compare with 5 egg scores
m4 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScoreFive, data=data[data$Generation=="G2",])
## Compare with EggScoreSmall
m5 <- lm(asin(sqrt(Emergence_Rate)) ~ Population + Test_environment + SA + Original_environment:Test_environment+EggScoreSmall, data=data[data$Generation=="G2",])
summary(m4)##
## Call:
## lm(formula = asin(sqrt(Emergence_Rate)) ~ Population + Test_environment +
## SA + Original_environment:Test_environment + EggScoreFive,
## data = data[data$Generation == "G2", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42998 -0.09038 -0.00473 0.07960 0.93896
##
## Coefficients: (3 not defined because of singularities)
## Estimate Std. Error
## (Intercept) 0.601300 0.055289
## PopulationBlackberry32 -0.377346 0.037959
## PopulationBlackberry33 -0.055881 0.036355
## PopulationBlackberry34 -0.228972 0.053250
## PopulationBlackberry35 -0.209601 0.036206
## PopulationBlackberry36 -0.107051 0.046078
## PopulationBlackberry37 -0.119387 0.038724
## PopulationBlackberry38 -0.320919 0.072830
## PopulationBlackberry39 -0.107259 0.041814
## PopulationBlackberry40 0.012021 0.039828
## PopulationBlackberry43 -0.022641 0.055390
## PopulationBlackberry44 -0.026557 0.039715
## PopulationBlackberry45 -0.236677 0.079817
## PopulationCherry103 0.044764 0.049636
## PopulationCherry104 -0.240816 0.052751
## PopulationCherry3 -0.481162 0.061106
## PopulationCherry47 -0.241190 0.042944
## PopulationCherry50 0.159812 0.057307
## PopulationCherry51 0.278492 0.155026
## PopulationCherry52 0.135919 0.083205
## PopulationCherry6 -0.360564 0.068615
## PopulationCherry7 -0.190717 0.043350
## PopulationStrawberry42 -0.007520 0.041778
## PopulationStrawberry44 -0.255064 0.044849
## PopulationStrawberry53 -0.136571 0.040989
## Test_environmentCherry 0.001734 0.047496
## Test_environmentStrawberry -0.117162 0.037091
## SA1 0.069558 0.037314
## EggScoreFive2 -0.016379 0.028899
## EggScoreFive3 -0.036700 0.030340
## EggScoreFive4 -0.086037 0.033727
## EggScoreFive5 -0.099671 0.047246
## Test_environmentBlackberry:Original_environmentCherry 0.003612 0.048933
## Test_environmentCherry:Original_environmentCherry -0.043605 0.059561
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.029804 0.064476
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## t value Pr(>|t|)
## (Intercept) 10.876 < 2e-16 ***
## PopulationBlackberry32 -9.941 < 2e-16 ***
## PopulationBlackberry33 -1.537 0.124831
## PopulationBlackberry34 -4.300 2.01e-05 ***
## PopulationBlackberry35 -5.789 1.17e-08 ***
## PopulationBlackberry36 -2.323 0.020518 *
## PopulationBlackberry37 -3.083 0.002148 **
## PopulationBlackberry38 -4.406 1.26e-05 ***
## PopulationBlackberry39 -2.565 0.010569 *
## PopulationBlackberry40 0.302 0.762906
## PopulationBlackberry43 -0.409 0.682871
## PopulationBlackberry44 -0.669 0.503964
## PopulationBlackberry45 -2.965 0.003151 **
## PopulationCherry103 0.902 0.367523
## PopulationCherry104 -4.565 6.13e-06 ***
## PopulationCherry3 -7.874 1.75e-14 ***
## PopulationCherry47 -5.616 3.06e-08 ***
## PopulationCherry50 2.789 0.005470 **
## PopulationCherry51 1.796 0.072960 .
## PopulationCherry52 1.634 0.102911
## PopulationCherry6 -5.255 2.10e-07 ***
## PopulationCherry7 -4.399 1.30e-05 ***
## PopulationStrawberry42 -0.180 0.857225
## PopulationStrawberry44 -5.687 2.07e-08 ***
## PopulationStrawberry53 -3.332 0.000919 ***
## Test_environmentCherry 0.037 0.970893
## Test_environmentStrawberry -3.159 0.001669 **
## SA1 1.864 0.062819 .
## EggScoreFive2 -0.567 0.571108
## EggScoreFive3 -1.210 0.226929
## EggScoreFive4 -2.551 0.011004 *
## EggScoreFive5 -2.110 0.035331 *
## Test_environmentBlackberry:Original_environmentCherry 0.074 0.941190
## Test_environmentCherry:Original_environmentCherry -0.732 0.464411
## Test_environmentStrawberry:Original_environmentCherry NA NA
## Test_environmentBlackberry:Original_environmentStrawberry 0.462 0.644079
## Test_environmentCherry:Original_environmentStrawberry NA NA
## Test_environmentStrawberry:Original_environmentStrawberry NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1484 on 567 degrees of freedom
## Multiple R-squared: 0.5051, Adjusted R-squared: 0.4754
## F-statistic: 17.02 on 34 and 567 DF, p-value: < 2.2e-16model.sel(m1, m2, m3, m4, m5)## Model selection table
## (Int) Ppl SA Tst_env Org_env:Tst_env log(Nb_egg+1) EgS ESF ESS df logLik
## m2 0.9142 + + + + -0.07308 33 318.394
## m3 0.6020 + + + + + 35 312.255
## m5 0.5991 + + + + + 37 314.175
## m4 0.6013 + + + + + 36 312.318
## m1 0.5794 + + + + 32 305.067
## AICc delta weight
## m2 -566.8 0.00 1
## m3 -550.1 16.78 0
## m5 -549.4 17.47 0
## m4 -547.9 18.92 0
## m1 -542.4 24.42 0
## Models ranked by AICc(x)## Cl= model m5 provides a better description of the data than model m4 (variation in the number of eggs between 50 and 100 is more important than between 150 and 250)2.5.6 Power Analysis
m2 <- lm(asin(sqrt(Emergence_Rate)) ~ -1 + Population + Test_environment + Original_environment:Test_environment, data=data[data$Generation=="G2",])
simdata <- function(seed=1, nrep = 2, npop = 3, varres=0){
#print(seed)
set.seed(seed)
## Sample populations
PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
## Create dataset
newdata <- expand.grid(Pop_new_name=1:(3*npop), Test_environment=levels(data$Test_environment), Rep=1:nrep)
newdata$Original_environment <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
## Add original names to get the right coef in the model
newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
newdata$Population <- newdata$Pop_new_name
levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
## Get design matrix
mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)
## Check that the order of the coefficients is the same
data.frame(colnames(mat), names(coef(m2))[names(coef(m2))%in%colnames(mat)])
coefsim <- coef(m2)[names(coef(m2))%in%colnames(mat)]
coefsim <- ifelse(is.na(coefsim), 0, coefsim)
newdata$resp <- mat%*%coefsim + rnorm(n=nrow(newdata), mean = 0, sd = varres)
m3 <- lm(resp ~ -1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)
data.frame(coefsim, coef(m3))
newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
m4 <- lm(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment, data=newdata)
summary(m4)
anova(m4)
## F test for SA
(Fratio_int = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
return(pvalue_int)
}
## With 6 populations per fruit
simdata(seed=2, nrep = 30, npop = 6)
simpval <- sapply(1:500, simdata, nrep = 10, npop = 6, varres=sigma(m2))
mean(simpval)
hist(simpval)
## Compute power
mean(simpval<0.05)
simpval10rep <- sapply(1:500, simdata, nrep = 10, npop = 3, varres=sigma(m2))
simpval20rep <- sapply(1:500, simdata, nrep = 20, npop = 3, varres=sigma(m2))
simpval25rep <- sapply(1:500, simdata, nrep = 25, npop = 3, varres=sigma(m2))
simpval30rep <- sapply(1:500, simdata, nrep = 30, npop = 3, varres=sigma(m2))
simpval35rep <- sapply(1:500, simdata, nrep = 35, npop = 3, varres=sigma(m2))
simpval40rep <- sapply(1:500, simdata, nrep = 40, npop = 3, varres=sigma(m2))
simpval60rep <- sapply(1:500, simdata, nrep = 60, npop = 3, varres=sigma(m2))
simpval80rep <- sapply(1:500, simdata, nrep = 80, npop = 3, varres=sigma(m2))
simpval100rep <- sapply(1:500, simdata, nrep = 100, npop = 3, varres=sigma(m2))
## Compute power
val <- data.frame(simpval10rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval60rep, simpval80rep, simpval100rep)
computepower <- function(x){
return(mean(x<0.05))
}
statpower <- apply(val, 2, computepower)
plot(c(10, 20, 25, 30, 35, 40, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power")3 Preference
3.1 Load data
setwd("/Users/rodenico/Documents/Pro/Articles/2021_NGS_EvolExp/Simulations")
data <- read.table("DATACOMPLET_PREF_POPNAT2018.csv", header=TRUE, sep=";")
head(data)## Generation Experiment Box Date_P Original_environment Pop_code Pop_name
## 1 G0 Plasticity 860 19/09/2018 Blackberry 33 Cimetiere
## 2 G0 Plasticity 860 19/09/2018 Blackberry 33 Cimetiere
## 3 G0 Plasticity 860 19/09/2018 Blackberry 33 Cimetiere
## 4 G0 Plasticity 860 19/09/2018 Blackberry 33 Cimetiere
## 5 G0 Plasticity 860 19/09/2018 Blackberry 33 Cimetiere
## 6 G0 Plasticity 860 19/09/2018 Blackberry 33 Cimetiere
## Pop_original_name Population Row Column Test_environment Nb_eggs Date_C_O
## 1 Blackberry33 Blackberry33 1 1 Cranberry 0 13/12/2018
## 2 Blackberry33 Blackberry33 1 2 Fig 0 13/12/2018
## 3 Blackberry33 Blackberry33 1 3 Raspberry 0 13/12/2018
## 4 Blackberry33 Blackberry33 1 4 Rosehips 0 13/12/2018
## 5 Blackberry33 Blackberry33 2 1 Kiwi 0 13/12/2018
## 6 Blackberry33 Blackberry33 2 2 Strawberry 1 13/12/2018
## Obs_O
## 1 CD
## 2 CD
## 3 CD
## 4 CD
## 5 CD
## 6 CDdata <- data[data$Generation=="G2",]
#Remove GF
data <- data[data$Test_environment!="GF",]
data$Test_environment <- factor(data$Test_environment)
## Remove WT3
data <- data[data$Population!="WT3",]
data$Original_environment <- factor(data$Original_environment)
## Remove Cherry52 which has a very high mean
#data <- data[data$Population!="Cherry52",]
data$Population <- factor(data$Population)
data$Nb_eggs <- as.numeric(as.character(data$Nb_eggs))
data$Row_Col <- as.factor(paste(data$Row, data$Column, sep="_"))
head(data)## Generation Experiment Box Date_P Original_environment Pop_code
## 2665 G2 Adaptation 999 14/11/2018 Blackberry 45
## 2666 G2 Adaptation 999 14/11/2018 Blackberry 45
## 2667 G2 Adaptation 999 14/11/2018 Blackberry 45
## 2668 G2 Adaptation 999 14/11/2018 Blackberry 45
## 2669 G2 Adaptation 999 14/11/2018 Blackberry 45
## 2670 G2 Adaptation 999 14/11/2018 Blackberry 45
## Pop_name Pop_original_name Population Row Column Test_environment
## 2665 Uzes Blackberry45 Blackberry45 1 1 Grape
## 2666 Uzes Blackberry45 Blackberry45 1 2 Blackberry
## 2667 Uzes Blackberry45 Blackberry45 1 3 Cherry
## 2668 Uzes Blackberry45 Blackberry45 1 4 Strawberry
## 2669 Uzes Blackberry45 Blackberry45 2 1 Fig
## 2670 Uzes Blackberry45 Blackberry45 2 2 Rosehips
## Nb_eggs Date_C_O Obs_O Row_Col
## 2665 0 29/11/2018 LO 1_1
## 2666 0 29/11/2018 LO 1_2
## 2667 16 29/11/2018 LO 1_3
## 2668 12 29/11/2018 LO 1_4
## 2669 0 29/11/2018 LO 2_1
## 2670 0 29/11/2018 LO 2_2sort(unique(data$Nb_eggs))## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## [19] 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
## [37] 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53
## [55] 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 70 71 72
## [73] 74 75 76 77 78 79 80 81 82 83 84 85 87 89 90 95 98 100
## [91] 101 107 110 111 115 124 125 128 130 139 141 1593.2 3 fruits
3.2.1 Explore data with blue
m0 <- lm(Nb_eggs~Original_environment*Test_environment, data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])
summary(m0)##
## Call:
## lm(formula = Nb_eggs ~ Original_environment * Test_environment,
## data = data[data$Test_environment == "Blackberry" | data$Test_environment ==
## "Cherry" | data$Test_environment == "Strawberry", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -30.385 -13.622 -5.925 8.260 134.075
##
## Coefficients:
## Estimate Std. Error
## (Intercept) 24.622 2.196
## Original_environmentCherry -8.526 3.730
## Original_environmentStrawberry -2.697 4.079
## Test_environmentCherry -1.204 3.106
## Test_environmentStrawberry -17.684 3.106
## Original_environmentCherry:Test_environmentCherry 15.493 5.275
## Original_environmentStrawberry:Test_environmentCherry 4.204 5.769
## Original_environmentCherry:Test_environmentStrawberry 12.261 5.275
## Original_environmentStrawberry:Test_environmentStrawberry 8.534 5.769
## t value Pr(>|t|)
## (Intercept) 11.212 < 2e-16 ***
## Original_environmentCherry -2.286 0.02263 *
## Original_environmentStrawberry -0.661 0.50871
## Test_environmentCherry -0.388 0.69840
## Test_environmentStrawberry -5.694 2.01e-08 ***
## Original_environmentCherry:Test_environmentCherry 2.937 0.00345 **
## Original_environmentStrawberry:Test_environmentCherry 0.729 0.46645
## Original_environmentCherry:Test_environmentStrawberry 2.324 0.02047 *
## Original_environmentStrawberry:Test_environmentStrawberry 1.479 0.13963
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 21.74 on 561 degrees of freedom
## Multiple R-squared: 0.1117, Adjusted R-squared: 0.09901
## F-statistic: 8.816 on 8 and 561 DF, p-value: 2.139e-11## Compute proportions
sumeggsperbox <- tapply(data$Nb_eggs, data$Box, sum)
data$sumeggsperbox <- sumeggsperbox[match(as.character(data$Box), names(sumeggsperbox))]
data$prop <- data$Nb_eggs/data$sumeggsperbox
## Check for local adaptation pattern in proportion
tapply(data$prop[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], list(data$Original_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], data$Test_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"]), mean)## Apricot Blackberry Blackcurrant Cherry Cranberry Fig Grape Kiwi
## Blackberry NA 0.12654832 NA 0.1309270 NA NA NA NA
## Cherry NA 0.08179375 NA 0.1952284 NA NA NA NA
## Strawberry NA 0.12029602 NA 0.1280855 NA NA NA NA
## Raspberry Rosehips Strawberry Tomato
## Blackberry NA NA 0.03872770 NA
## Cherry NA NA 0.06151998 NA
## Strawberry NA NA 0.08369925 NA## Check for local adaptation pattern in te nulber of eggs
## By pop
tapply(data$Nb_eggs[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], list(data$Population[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], data$Test_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"]), mean)## Apricot Blackberry Blackcurrant Cherry Cranberry Fig Grape Kiwi
## Blackberry31 NA 31.250000 NA 24.250000 NA NA NA NA
## Blackberry32 NA 12.000000 NA 13.307692 NA NA NA NA
## Blackberry33 NA 28.000000 NA 32.583333 NA NA NA NA
## Blackberry34 NA 10.000000 NA 6.333333 NA NA NA NA
## Blackberry35 NA 34.857143 NA 23.571429 NA NA NA NA
## Blackberry36 NA 37.800000 NA 35.200000 NA NA NA NA
## Blackberry37 NA 25.181818 NA 22.818182 NA NA NA NA
## Blackberry38 NA 55.000000 NA 37.000000 NA NA NA NA
## Blackberry39 NA 12.000000 NA 29.714286 NA NA NA NA
## Blackberry40 NA 15.714286 NA 9.000000 NA NA NA NA
## Blackberry43 NA 9.333333 NA 12.000000 NA NA NA NA
## Blackberry44 NA 31.538462 NA 30.846154 NA NA NA NA
## Blackberry45 NA 0.000000 NA 16.000000 NA NA NA NA
## Cherry103 NA 11.000000 NA 23.142857 NA NA NA NA
## Cherry104 NA 8.800000 NA 17.000000 NA NA NA NA
## Cherry3 NA 8.333333 NA 25.000000 NA NA NA NA
## Cherry47 NA 16.076923 NA 33.307692 NA NA NA NA
## Cherry50 NA 35.500000 NA 37.250000 NA NA NA NA
## Cherry52 NA 0.500000 NA 28.500000 NA NA NA NA
## Cherry6 NA 53.250000 NA 57.250000 NA NA NA NA
## Cherry7 NA 9.000000 NA 27.857143 NA NA NA NA
## Strawberry42 NA 21.125000 NA 54.250000 NA NA NA NA
## Strawberry44 NA 16.444444 NA 25.333333 NA NA NA NA
## Strawberry53 NA 24.347826 NA 14.565217 NA NA NA NA
## Raspberry Rosehips Strawberry Tomato
## Blackberry31 NA NA 12.2500000 NA
## Blackberry32 NA NA 2.0000000 NA
## Blackberry33 NA NA 8.1666667 NA
## Blackberry34 NA NA 5.3333333 NA
## Blackberry35 NA NA 4.3571429 NA
## Blackberry36 NA NA 15.2000000 NA
## Blackberry37 NA NA 5.0909091 NA
## Blackberry38 NA NA 5.0000000 NA
## Blackberry39 NA NA 4.4285714 NA
## Blackberry40 NA NA 9.5714286 NA
## Blackberry43 NA NA 5.0000000 NA
## Blackberry44 NA NA 9.1538462 NA
## Blackberry45 NA NA 12.0000000 NA
## Cherry103 NA NA 0.4285714 NA
## Cherry104 NA NA 3.6000000 NA
## Cherry3 NA NA 28.6666667 NA
## Cherry47 NA NA 8.1538462 NA
## Cherry50 NA NA 17.5000000 NA
## Cherry52 NA NA 1.5000000 NA
## Cherry6 NA NA 40.2500000 NA
## Cherry7 NA NA 7.7142857 NA
## Strawberry42 NA NA 20.8750000 NA
## Strawberry44 NA NA 8.4444444 NA
## Strawberry53 NA NA 11.6521739 NA ## By original environment
tapply(data$Nb_eggs[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], list(data$Original_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], data$Test_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"]), mean)## Apricot Blackberry Blackcurrant Cherry Cranberry Fig Grape Kiwi
## Blackberry NA 24.62245 NA 23.41837 NA NA NA NA
## Cherry NA 16.09615 NA 30.38462 NA NA NA NA
## Strawberry NA 21.92500 NA 24.92500 NA NA NA NA
## Raspberry Rosehips Strawberry Tomato
## Blackberry NA NA 6.938776 NA
## Cherry NA NA 10.673077 NA
## Strawberry NA NA 12.775000 NAtapply(data$Nb_eggs[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], list(data$Original_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"], data$Test_environment[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"]), var)## Apricot Blackberry Blackcurrant Cherry Cranberry Fig Grape Kiwi
## Blackberry NA 657.9694 NA 515.5036 NA NA NA NA
## Cherry NA 381.8533 NA 623.6531 NA NA NA NA
## Strawberry NA 722.4814 NA 836.4814 NA NA NA NA
## Raspberry Rosehips Strawberry Tomato
## Blackberry NA NA 85.29518 NA
## Cherry NA NA 405.36161 NA
## Strawberry NA NA 264.28141 NArange(data$Nb_eggs)## [1] 0 159## Check for the presence of negative correlations after accounting for the heterogeneity between populations (i.e. some lay few or a lot of eggs)
m1 <- lm(log(Nb_eggs+1) ~ Population, data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])
data$resid <- NA
data$resid[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry"] <- residuals(m1)
meanbypopbytestenv <- as.data.frame(tapply(data$resid[data$Generation=="G2"], list(data$Population[data$Generation=="G2"], data$Test_environment[data$Generation=="G2"]), mean))
## Cherry ~ Blackberry
plot(meanbypopbytestenv$Blackberry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Cherry", pch=16)
legend("bottomleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Strawberry ~ Cherry
plot(meanbypopbytestenv$Cherry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Cherry", ylab="Strawberry", pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Strawberry ~ Blackberry
plot(meanbypopbytestenv$Blackberry~ meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), xlab="Blackberry", ylab="Strawberry", pch=16)
legend("bottomleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
3.2.2 Explore data with BLUPS
m0 <- lmer(log(Nb_eggs+1)~ Population + Test_environment + (0 + Test_environment | Population), data = data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])
## Correlation between fixed effects
plot(coef(m1), fixef(m0))
summary(m0)
BLUP <- ranef(m0)$Population
pairs(BLUP)
## Positive correlation between blups and blues, as expected
plot(BLUP$Test_environmentBlackberry, meanbypopbytestenv$Blackberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Positive correlation betwen blups and blues, as expected
plot(BLUP$Test_environmentCherry, meanbypopbytestenv$Cherry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("topleft", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)
## Negative correlation betwen blups and blues, weird
plot(BLUP$Test_environmentStrawberry, meanbypopbytestenv$Strawberry, col=as.numeric(data$Original_environment[match(rownames(meanbypopbytestenv), data$Population)]), pch=16)
legend("bottomright", levels(data$Original_environment), col=as.numeric(as.factor(levels(data$Original_environment))), pch=16)3.2.3 Check that model with box random effect is OK for the analyses for analyses
data$SA <- as.factor(ifelse(as.character(data$Original_environment)==as.character(data$Test_environment), 1, 0))
data$Box <- as.factor(data$Box)
data$Obs <- as.factor(1:nrow(data))
npop=3
PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
m1 <- glmer(Nb_eggs ~ 0 + Original_environment + Test_environment + SA + Original_environment:Test_environment+ (1|Box)+(1|Obs), data=data[(data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry")&data$Population%in%c(PopBlackberry, PopCherry, PopStrawberry),], family="poisson")## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00691162 (tol = 0.002, component 1)summary(m1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Nb_eggs ~ 0 + Original_environment + Test_environment + SA +
## Original_environment:Test_environment + (1 | Box) + (1 | Obs)
## Data:
## data[(data$Test_environment == "Blackberry" | data$Test_environment ==
## "Cherry" | data$Test_environment == "Strawberry") & data$Population %in%
## c(PopBlackberry, PopCherry, PopStrawberry), ]
##
## AIC BIC logLik deviance df.resid
## 1110.4 1143.3 -544.2 1088.4 136
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.11291 -0.38074 0.02317 0.10502 0.46559
##
## Random effects:
## Groups Name Variance Std.Dev.
## Obs (Intercept) 1.5308 1.2372
## Box (Intercept) 0.5562 0.7458
## Number of obs: 147, groups: Obs, 147; Box, 49
##
## Fixed effects:
## Estimate Std. Error
## Original_environmentBlackberry 1.14738 0.36282
## Original_environmentCherry 1.65103 0.48325
## Original_environmentStrawberry 2.34639 0.36725
## Test_environmentCherry 0.56461 0.43532
## Test_environmentStrawberry -0.73526 0.31725
## SA1 0.52175 0.31750
## Original_environmentCherry:Test_environmentCherry -0.08221 0.85490
## Original_environmentStrawberry:Test_environmentCherry 0.45414 0.53071
## Original_environmentCherry:Test_environmentStrawberry 0.24455 0.67875
## z value Pr(>|z|)
## Original_environmentBlackberry 3.162 0.001565 **
## Original_environmentCherry 3.417 0.000634 ***
## Original_environmentStrawberry 6.389 1.67e-10 ***
## Test_environmentCherry 1.297 0.194625
## Test_environmentStrawberry -2.318 0.020473 *
## SA1 1.643 0.100320
## Original_environmentCherry:Test_environmentCherry -0.096 0.923393
## Original_environmentStrawberry:Test_environmentCherry 0.856 0.392154
## Original_environmentCherry:Test_environmentStrawberry 0.360 0.718623
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## Orgn_B Orgn_C Orgn_S Tst_nC Tst_nS SA1 O_C:T_C O_S:T_
## Orgnl_nvrnC 0.023
## Orgnl_nvrnS 0.393 0.010
## Tst_nvrnmnC -0.650 -0.008 -0.321
## Tst_nvrnmnS -0.393 0.002 -0.433 0.329
## SA1 -0.491 -0.006 -0.437 0.405 0.005
## Orgnl_C:T_C 0.504 -0.432 0.322 -0.657 -0.170 -0.576
## Orgnl_S:T_C 0.262 0.000 -0.261 -0.598 0.030 -0.030 0.316
## Orgnl_C:T_S 0.180 -0.547 0.201 -0.152 -0.468 -0.001 0.385 -0.014
## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00691162 (tol = 0.002, component 1)m1bis <- glmer(Nb_eggs ~ 0 + Original_environment + Test_environment + (1|Box)+(1|Obs), data=data[(data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry")&data$Population%in%c(PopBlackberry, PopCherry, PopStrawberry),], family="poisson")
summary(m1bis)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: Nb_eggs ~ 0 + Original_environment + Test_environment + (1 |
## Box) + (1 | Obs)
## Data:
## data[(data$Test_environment == "Blackberry" | data$Test_environment ==
## "Cherry" | data$Test_environment == "Strawberry") & data$Population %in%
## c(PopBlackberry, PopCherry, PopStrawberry), ]
##
## AIC BIC logLik deviance df.resid
## 1106.4 1127.4 -546.2 1092.4 140
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.14374 -0.37310 0.02474 0.10761 0.40015
##
## Random effects:
## Groups Name Variance Std.Dev.
## Obs (Intercept) 1.6136 1.2703
## Box (Intercept) 0.5291 0.7274
## Number of obs: 147, groups: Obs, 147; Box, 49
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## Original_environmentBlackberry 1.2851 0.2959 4.343 1.41e-05 ***
## Original_environmentCherry 1.8466 0.3710 4.978 6.43e-07 ***
## Original_environmentStrawberry 2.6376 0.3051 8.644 < 2e-16 ***
## Test_environmentCherry 0.6224 0.2768 2.248 0.0246 *
## Test_environmentStrawberry -0.6923 0.2866 -2.415 0.0157 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## Orgn_B Orgn_C Orgn_S Tst_nC
## Orgnl_nvrnC 0.269
## Orgnl_nvrnS 0.305 0.251
## Tst_nvrnmnC -0.496 -0.406 -0.475
## Tst_nvrnmnS -0.451 -0.377 -0.462 0.497simdata <- function(seed=1, nrep = 1000, npop = 3, SA=FALSE){
set.seed(seed)
## Create dataset
newdata <- expand.grid(Original_environment=levels(data$Original_environment), Test_environment=c("Blackberry", "Cherry", "Strawberry"), Box=1:nrep)
newdata$SA <- as.factor(ifelse(as.character(newdata$Original_environment)==as.character(newdata$Test_environment), 1, 0))
newdata$Box <- as.factor(newdata$Box)
newdata$Obs <- as.factor(1:nrow(newdata))
if(SA==TRUE){
## Get design matrix
mat <- model.matrix(~0 + Original_environment + Test_environment + SA + Original_environment:Test_environment , data=newdata)
m <- m1
## Drop Strawberry x Strawberry interaction which is not identifiable
mat <- mat[, -ncol(mat)]
}else{
mat <- model.matrix(~0 + Original_environment + Test_environment , data=newdata)
m <- m1bis
}
## Check that the order of the coefficients is the same
data.frame(colnames(mat), names(fixef(m))[names(fixef(m))%in%colnames(mat)])
## Get the coefficient of the right populations
coefsim <- fixef(m)[names(fixef(m))%in%colnames(mat)]
coefsim <- ifelse(is.na(coefsim), 0, coefsim)
## Add box random effect
mat2 <- model.matrix(~0 + Box, data=newdata)
coefrand <- rnorm(n=nrep, sd=sqrt(data.frame(VarCorr(m))$vcov[2]))
## Simulate data
newdata$resp <- exp(mat%*%coefsim + mat2%*%coefrand + rnorm(n=nrow(newdata), mean = 0, sd = sqrt(as.data.frame(VarCorr(m))$vcov[2])))
newdata$Nb_eggs <- sapply(newdata$resp, rpois, n=1)
if(SA==TRUE){
## Check fit of interaction model
m3 <- glmer(Nb_eggs ~ 0 + Original_environment + Test_environment + SA + Original_environment:Test_environment + (1|Box) + (1|Obs), data=newdata, family="poisson")
data.frame(coefsim, fixef(m3))
}
m4 <- lmer(log(Nb_eggs+1) ~ 0 + Original_environment + Test_environment + SA + Original_environment:Test_environment + (1|Box) , data=newdata)
anova(m4)
## F test for SA
(Fratio_mixed = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
## Compute P-value
(pvalue_mixed = 1 - pf(Fratio_mixed, 1, anova(m4)[4,1]))
## Fit SA model
m5 <- lm(log(Nb_eggs+1) ~ Original_environment + Test_environment + SA + Original_environment:Test_environment + Box, data=newdata)
summary(m5)
anova(m5)
## F test for SA
(Fratio_fixed = (anova(m5)[3,2]/anova(m5)[5,2])/(1/anova(m5)[5,1]))
## Compute P-value
(pvalue_fixed = 1 - pf(Fratio_fixed, 1, anova(m5)[5,1]))
return(c(Fratio_mixed=Fratio_mixed, pvalue_mixed=pvalue_mixed, Fratio_fixed=Fratio_fixed, pvalue_fixed=pvalue_fixed))
}
simpval <- sapply(1:100, simdata, nrep = 30, npop = 3, SA=TRUE)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.209911 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.030907 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.137669 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.181223 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0898904 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00268656 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.24508 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0842592 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0526515 (tol = 0.002, component 1)## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0970607 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00440259 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0135188 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0193863 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00596128 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.012778 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0459003 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00562554 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0624345 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0265195 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00780001 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0320825 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0383621 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00834037 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0663908 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0281588 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0544856 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0595906 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.048528 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.057708 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0494004 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0388851 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0379331 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0139703 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0831132 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0470658 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0554048 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0607758 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0552343 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0163608 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00481468 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0161518 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00475671 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00931366 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0624666 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0612396 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00323709 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.159664 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00256289 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0713899 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.134106 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00527572 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0245835 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.056269 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.10428 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0534692 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0784539 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0135546 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0571271 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.05802 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0536986 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0403313 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0586634 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.057353 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00429039 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0662628 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0454432 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0245136 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.121999 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0634159 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0576719 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0608196 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00231643 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00679617 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0377365 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0457307 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0580104 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0192054 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0220035 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0328747 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0576076 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0592644 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0615031 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0521504 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0543718 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.238321 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.105513 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.156114 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0157434 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0206585 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model failed to converge with max|grad| = 0.0582537 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0117851 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0335173 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0452748 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0288024 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00378826 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0661943 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0520898 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00386292 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0375417 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.0219405 (tol = 0.002, component 1)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## Power of 51%
mean(simpval[2,]<0.05)## [1] 0.12## Check false positive rate of 5%
simpval <- sapply(1:100, simdata, nrep = 100, npop = 3, SA=FALSE)## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient
## fixed-effect model matrix is rank deficient so dropping 1 column / coefficient## False positive rate of 5%
mean(simpval[2,]<0.05)## [1] 0.023.2.4 Analysis
data$SA <- as.factor(ifelse(as.character(data$Original_environment)==as.character(data$Test_environment), 1, 0))
data$Box <- as.factor(data$Box)
## Model with a box effect
m1 <- lmer(log(Nb_eggs+1) ~ Population + Test_environment + SA + (1|Box) + Original_environment:Test_environment, data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients## F test for SA with Box effect
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))## [1] 11.54878(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.04251289## Compute Rsquare
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))## [1] 0.7250628## Model with a box effect and a rwo/col effect
m2 <- lmer(log(Nb_eggs+1) ~ Population + Test_environment + SA + (1|Box) + (1|Row_Col) + Original_environment:Test_environment, data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients## F test for SA with Box effect
(Fratio_int = (anova(m2)[3,2]/anova(m2)[4,2])/(1/anova(m2)[4,1]))## [1] 12.46661(pvalue_int = 1 - pf(Fratio_int, 1, anova(m2)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.03861686## Compute Rsquare
(rsq2 <- 1-anova(m2)[4, 3]/((anova(m2)[3, 2]+anova(m2)[4, 2])/(anova(m2)[3, 1]+anova(m2)[4, 1])))## [1] 0.74137833.2.5 Power Analysis
m2 <- lmer(log(Nb_eggs+1) ~ -1 + Population + Test_environment + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=data[data$Test_environment=="Blackberry"|data$Test_environment=="Cherry"|data$Test_environment=="Strawberry",])
fixef(m2)
simdata <- function(seed=1, nrep = 2, npop = 3){
#print(seed)
set.seed(seed)
## Sample populations
PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
## Create dataset
newdata <- expand.grid(Box=factor(1:nrep), Pop_new_name=1:(3*npop), Test_environment=rep(c("Blackberry", "Cherry", "Strawberry"), each=4))
newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
newdata$Original_environment <- newdata$Pop_new_name
levels(newdata$Original_environment) <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
newdata$Box <- as.factor(paste(newdata$Box, newdata$Pop_new_name, sep="_"))
#newdata[newdata$Box==2&newdata$Pop_new_name==1,]
#newdata[newdata$Box=="1_1",]
#newdata <- newdata[order(newdata$Box),]
#newdata[newdata$Original_environment=="Blackberry",]
nlevels(newdata$Box)
## Add original names to get the right coef in the model
newdata$Population <- newdata$Pop_new_name
levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
unique(newdata$Population)
newdata$Row_Col <- NA
for (i in levels(newdata$Box)){
newdata$Row_Col[newdata$Box==i] <- sample(x=1:12, size=12, replace = FALSE)
}
newdata$Row_Col <- as.factor(newdata$Row_Col)
## Get design matrix
mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)
## Remove levels that are not identifiable
mat <- mat[, !colnames(mat)%in%c("Test_environmentStrawberry:Original_environmentCherry", "Test_environmentStrawberry:Original_environmentStrawberry")]
## Check that the order of the coefficients is the same
data.frame(colnames(mat), names(fixef(m2))[names(fixef(m2))%in%colnames(mat)])
## Get the coefficient of the right populations
coefsim <- fixef(m2)[names(fixef(m2))%in%colnames(mat)]
#coefsim <- ifelse(is.na(coefsim), 0, coefsim)
## Add box random effect
mat2 <- model.matrix(~0 + Box, data=newdata)
coefrand <- rnorm(n=nrep*npop*3, sd=sqrt(data.frame(VarCorr(m2))$vcov[1]))
## Add row/col random effect
mat3 <- model.matrix(~0 + Row_Col, data=newdata)
coefrand2 <- rnorm(n=12, sd=sqrt(data.frame(VarCorr(m2))$vcov[2]))
## Simulate data
newdata$resp <- mat%*%coefsim + mat2%*%coefrand + mat3%*%coefrand2 + rnorm(n=nrow(newdata), mean = 0, sd = sqrt(as.data.frame(VarCorr(m2))$vcov[2]))
## Fit interaction model
m3 <- lmer(resp ~ -1 + Population + Test_environment + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=newdata)
data.frame(coefsim, fixef(m3))
## Fit SA model assuming 9 different populations
newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
m4 <- lmer(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=newdata)
summary(m4)
anova(m4)
## F test for SA
(Fratio_int = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
## Compute P-value
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[4,1]))
return(pvalue_int)
}
simdata(seed=2, nrep = 30, npop = 3)
simpval <- sapply(1:500, simdata, nrep = 10, npop = 3)
mean(simpval)
## Compute power
mean(simpval<0.05)
hist(simpval)
nsim <- 500
simpval10rep <- sapply(1:nsim, simdata, nrep = 10, npop = 3)
simpval15rep <- sapply(1:nsim, simdata, nrep = 15, npop = 3)
simpval20rep <- sapply(1:nsim, simdata, nrep = 20, npop = 3)
simpval25rep <- sapply(1:nsim, simdata, nrep = 25, npop = 3)
simpval30rep <- sapply(1:nsim, simdata, nrep = 30, npop = 3)
simpval35rep <- sapply(1:nsim, simdata, nrep = 35, npop = 3)
simpval40rep <- sapply(1:nsim, simdata, nrep = 40, npop = 3)
simpval60rep <- sapply(1:nsim, simdata, nrep = 60, npop = 3)
simpval80rep <- sapply(1:nsim, simdata, nrep = 80, npop = 3)
simpval100rep <- sapply(1:nsim, simdata, nrep = 100, npop = 3)
## Compute power
val <- data.frame(simpval10rep, simpval15rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval60rep, simpval80rep, simpval100rep)
computepower <- function(x){
return(mean(x<0.05))
}
statpower <- apply(val, 2, computepower)
plot(c(10, 15, 20, 25, 30, 35, 40, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power", main="Preference (3 fruits)")3.3 12 fruits
3.3.1 Analysis
data$SA <- as.factor(ifelse(as.character(data$Original_environment)==as.character(data$Test_environment), 1, 0))
m1 <- lmer(log(Nb_eggs+1) ~ Population + Test_environment + SA + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=data[data$Generation=="G2",])## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficientssummary(m1)## Linear mixed model fit by REML ['lmerMod']
## Formula:
## log(Nb_eggs + 1) ~ Population + Test_environment + SA + Original_environment:Test_environment +
## (1 | Box) + (1 | Row_Col)
## Data: data[data$Generation == "G2", ]
##
## REML criterion at convergence: 7194.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.74710 -0.70531 0.00879 0.68438 2.89632
##
## Random effects:
## Groups Name Variance Std.Dev.
## Box (Intercept) 0.38239 0.6184
## Row_Col (Intercept) 0.03124 0.1768
## Residual 1.18303 1.0877
## Number of obs: 2280, groups: Box, 190; Row_Col, 12
##
## Fixed effects:
## Estimate
## (Intercept) 1.506391
## PopulationBlackberry32 -0.702176
## PopulationBlackberry33 0.391290
## PopulationBlackberry34 -0.538466
## PopulationBlackberry35 -0.065216
## PopulationBlackberry36 0.727676
## PopulationBlackberry37 -0.172517
## PopulationBlackberry38 0.188417
## PopulationBlackberry39 -0.152217
## PopulationBlackberry40 -0.407543
## PopulationBlackberry43 0.143136
## PopulationBlackberry44 0.415712
## PopulationBlackberry45 -1.315642
## PopulationCherry103 -0.834328
## PopulationCherry104 -0.820710
## PopulationCherry3 0.213004
## PopulationCherry47 -0.286906
## PopulationCherry50 -0.066922
## PopulationCherry52 -0.275113
## PopulationCherry6 0.951440
## PopulationCherry7 -0.364385
## PopulationStrawberry42 0.841909
## PopulationStrawberry44 -0.164577
## PopulationStrawberry53 -0.345964
## Test_environmentBlackberry 0.525783
## Test_environmentBlackcurrant 0.761040
## Test_environmentCherry 1.237440
## Test_environmentCranberry -0.065099
## Test_environmentFig 0.479561
## Test_environmentGrape 1.149459
## Test_environmentKiwi 0.726754
## Test_environmentRaspberry 0.671069
## Test_environmentRosehips 0.356897
## Test_environmentStrawberry 0.013150
## Test_environmentTomato 0.330173
## SA1 0.560535
## Test_environmentApricot:Original_environmentCherry 0.558490
## Test_environmentBlackberry:Original_environmentCherry 0.235338
## Test_environmentBlackcurrant:Original_environmentCherry 0.212224
## Test_environmentCherry:Original_environmentCherry -0.007308
## Test_environmentCranberry:Original_environmentCherry 0.415566
## Test_environmentFig:Original_environmentCherry 0.307003
## Test_environmentGrape:Original_environmentCherry -0.658191
## Test_environmentKiwi:Original_environmentCherry -0.146252
## Test_environmentRaspberry:Original_environmentCherry 0.130497
## Test_environmentRosehips:Original_environmentCherry 0.266882
## Test_environmentStrawberry:Original_environmentCherry 0.309021
## Test_environmentApricot:Original_environmentStrawberry 0.418854
## Test_environmentBlackberry:Original_environmentStrawberry 0.474579
## Test_environmentBlackcurrant:Original_environmentStrawberry 0.184944
## Test_environmentCherry:Original_environmentStrawberry -0.123624
## Test_environmentCranberry:Original_environmentStrawberry -0.006736
## Test_environmentFig:Original_environmentStrawberry -0.121108
## Test_environmentGrape:Original_environmentStrawberry -0.685244
## Test_environmentKiwi:Original_environmentStrawberry -0.409178
## Test_environmentRaspberry:Original_environmentStrawberry 0.156703
## Test_environmentRosehips:Original_environmentStrawberry 0.256809
## Std. Error t value
## (Intercept) 0.271768 5.543
## PopulationBlackberry32 0.311642 -2.253
## PopulationBlackberry33 0.316550 1.236
## PopulationBlackberry34 0.469519 -1.147
## PopulationBlackberry35 0.307372 -0.212
## PopulationBlackberry36 0.395370 1.840
## PopulationBlackberry37 0.322254 -0.535
## PopulationBlackberry38 0.735595 0.256
## PopulationBlackberry39 0.358934 -0.424
## PopulationBlackberry40 0.358934 -1.135
## PopulationBlackberry43 0.469519 0.305
## PopulationBlackberry44 0.311642 1.334
## PopulationBlackberry45 0.735595 -1.789
## PopulationCherry103 0.401009 -2.081
## PopulationCherry104 0.433927 -1.891
## PopulationCherry3 0.502417 0.424
## PopulationCherry47 0.359298 -0.799
## PopulationCherry50 0.460805 -0.145
## PopulationCherry52 0.576702 -0.477
## PopulationCherry6 0.460805 2.065
## PopulationCherry7 0.355601 -1.025
## PopulationStrawberry42 0.398242 2.114
## PopulationStrawberry44 0.389764 -0.422
## PopulationStrawberry53 0.345523 -1.001
## Test_environmentBlackberry 0.328650 1.600
## Test_environmentBlackcurrant 0.155951 4.880
## Test_environmentCherry 0.155700 7.948
## Test_environmentCranberry 0.155495 -0.419
## Test_environmentFig 0.155674 3.081
## Test_environmentGrape 0.155742 7.381
## Test_environmentKiwi 0.155801 4.665
## Test_environmentRaspberry 0.155733 4.309
## Test_environmentRosehips 0.156020 2.288
## Test_environmentStrawberry 0.156214 0.084
## Test_environmentTomato 0.155728 2.120
## SA1 0.289287 1.938
## Test_environmentApricot:Original_environmentCherry 0.264159 2.114
## Test_environmentBlackberry:Original_environmentCherry 0.421525 0.558
## Test_environmentBlackcurrant:Original_environmentCherry 0.264235 0.803
## Test_environmentCherry:Original_environmentCherry 0.359791 -0.020
## Test_environmentCranberry:Original_environmentCherry 0.264373 1.572
## Test_environmentFig:Original_environmentCherry 0.264138 1.162
## Test_environmentGrape:Original_environmentCherry 0.264271 -2.491
## Test_environmentKiwi:Original_environmentCherry 0.264306 -0.553
## Test_environmentRaspberry:Original_environmentCherry 0.264269 0.494
## Test_environmentRosehips:Original_environmentCherry 0.264409 1.009
## Test_environmentStrawberry:Original_environmentCherry 0.264158 1.170
## Test_environmentApricot:Original_environmentStrawberry 0.289341 1.448
## Test_environmentBlackberry:Original_environmentStrawberry 0.501087 0.947
## Test_environmentBlackcurrant:Original_environmentStrawberry 0.289304 0.639
## Test_environmentCherry:Original_environmentStrawberry 0.289450 -0.427
## Test_environmentCranberry:Original_environmentStrawberry 0.289217 -0.023
## Test_environmentFig:Original_environmentStrawberry 0.289228 -0.419
## Test_environmentGrape:Original_environmentStrawberry 0.288957 -2.371
## Test_environmentKiwi:Original_environmentStrawberry 0.289304 -1.414
## Test_environmentRaspberry:Original_environmentStrawberry 0.289345 0.542
## Test_environmentRosehips:Original_environmentStrawberry 0.289024 0.889##
## Correlation matrix not shown by default, as p = 57 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need it## fit warnings:
## fixed-effect model matrix is rank deficient so dropping 3 columns / coefficientsanova(m1)## Analysis of Variance Table
## npar Sum Sq Mean Sq F value
## Population 23 90.908 3.9525 3.3410
## Test_environment 11 290.703 26.4276 22.3388
## SA 1 18.790 18.7904 15.8832
## Test_environment:Original_environment 21 55.756 2.6551 2.2443## F test for SA with Box effect
(Fratio_int = (anova(m1)[3,2]/anova(m1)[4,2])/(1/anova(m1)[4,1]))## [1] 7.077208(pvalue_int = 1 - pf(Fratio_int, 1, anova(m1)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)## [1] 0.01464063## Compute Rsquare
(rsq2 <- 1-anova(m1)[4, 3]/((anova(m1)[3, 2]+anova(m1)[4, 2])/(anova(m1)[3, 1]+anova(m1)[4, 1])))## [1] 0.21644633.3.2 Power Analysis
m2 <- lmer(log(Nb_eggs+1) ~ -1+Population + Test_environment + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=data[data$Generation=="G2",])
fixef(m2)
hist(log(data$Nb_eggs[data$Generation=="G2"]+1))
simdata <- function(seed=1, nrep = 2, npop = 3){
set.seed(seed)
## Sample populations
PopBlackberry <- sort(sample(levels(data$Population)[grep("Blackberry", levels(data$Population))], size=npop, replace = TRUE))
PopCherry <- sort(sample(levels(data$Population)[grep("Cherry", levels(data$Population))], size=npop, replace = TRUE))
PopStrawberry <- sort(sample(levels(data$Population)[grep("Strawberry", levels(data$Population))], size=npop, replace = TRUE))
## Create dataset
newdata <- expand.grid(Pop_new_name=1:(3*npop), Test_environment=levels(data$Test_environment), Box=factor(1:nrep))
newdata$Original_environment <- rep(c("Blackberry", "Cherry", "Strawberry"), each=npop)
## Add original names to get the right coef in the model
newdata$Pop_new_name <- as.factor(newdata$Pop_new_name)
newdata$Population <- newdata$Pop_new_name
levels(newdata$Population) <- c(PopBlackberry, PopCherry, PopStrawberry)
newdata$Box <- as.factor(paste(newdata$Box, newdata$Pop_new_name, sep="_"))
#newdata[newdata$Pop_new_name=="1",]
#newdata[order(newdata$Box),]
#nlevels(newdata$Box)
unique(newdata$Population)
newdata$Row_Col <- NA
for (i in levels(newdata$Box)){
newdata$Row_Col[newdata$Box==i] <- sample(x=1:12, size=12, replace = FALSE)
}
newdata$Row_Col <- as.factor(newdata$Row_Col)
## Get design matrix
mat <- model.matrix(~-1 + Population + Test_environment + Original_environment:Test_environment, data=newdata)
## Remove levels that are not identifiable
mat <- mat[, !colnames(mat)%in%colnames(mat)[!colnames(mat)%in%names(fixef(m2))]]
## Check that the order of the coefficients is the same
data.frame(colnames(mat), names(fixef(m2))[names(fixef(m2))%in%colnames(mat)])
coefsim <- fixef(m2)[names(fixef(m2))%in%colnames(mat)]
#coefsim <- ifelse(is.na(coefsim), 0, coefsim)
## Add box random effect
mat2 <- model.matrix(~0 + Box, data=newdata)
coefrand <- rnorm(n=nrep*npop*3, sd=sqrt(data.frame(VarCorr(m2))$vcov[1]))
## Add row/col random effect
mat3 <- model.matrix(~0 + Row_Col, data=newdata)
coefrand2 <- rnorm(n=12, sd=sqrt(data.frame(VarCorr(m2))$vcov[2]))
## Simulate data
newdata$resp <- mat%*%coefsim + mat2%*%coefrand + mat3%*%coefrand2 + rnorm(n=nrow(newdata), mean = 0, sd = sqrt(as.data.frame(VarCorr(m2))$vcov[3]))
m3 <- lmer(resp ~ -1 + Population + Test_environment + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=newdata)
data.frame(coefsim, fixef(m3))
newdata$SA <- as.factor(ifelse(newdata$Original_environment==newdata$Test_environment, 1, 0))
m4 <- lmer(resp ~ Pop_new_name + Test_environment + SA + Original_environment:Test_environment + (1|Box) + (1|Row_Col), data=newdata)
summary(m4)
anova(m4)
## F test for SA
(Fratio_int = (anova(m4)[3,2]/anova(m4)[4,2])/(1/anova(m4)[4,1]))
(pvalue_int = 1 - pf(Fratio_int, 1, anova(m4)[4,1])) #the correct test (see equation D7 in Appendix D of the paper)
return(pvalue_int)
}
simdata(seed=2, nrep = 30, npop = 3)
simpval <- sapply(1:500, simdata, nrep = 10, npop = 3)
mean(simpval)
## Compute power
mean(simpval<0.05)
hist(simpval)
nsim <- 500
simpval10rep <- sapply(1:nsim, simdata, nrep = 10, npop = 3)
simpval15rep <- sapply(1:nsim, simdata, nrep = 15, npop = 3)
simpval20rep <- sapply(1:nsim, simdata, nrep = 20, npop = 3)
simpval25rep <- sapply(1:nsim, simdata, nrep = 25, npop = 3)
simpval30rep <- sapply(1:nsim, simdata, nrep = 30, npop = 3)
simpval35rep <- sapply(1:nsim, simdata, nrep = 35, npop = 3)
simpval40rep <- sapply(1:nsim, simdata, nrep = 40, npop = 3)
simpval60rep <- sapply(1:nsim, simdata, nrep = 60, npop = 3)
simpval80rep <- sapply(1:nsim, simdata, nrep = 80, npop = 3)
simpval100rep <- sapply(1:nsim, simdata, nrep = 100, npop = 3)
## Compute power
val <- data.frame(simpval10rep, simpval15rep, simpval20rep, simpval25rep, simpval30rep, simpval35rep, simpval40rep, simpval60rep, simpval80rep, simpval100rep)
computepower <- function(x){
return(mean(x<0.05))
}
statpower <- apply(val, 2, computepower)
plot(c(10, 15, 20, 25, 30, 35, 40, 60, 80, 100), statpower, las=1, xlab="Number of replicates", ylab="Statistical Power", main="Preference (12 fruits)")